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ABSTRACT 

The secular evolution of an infinitely thin tepid isolated galactic disc made of a finite number of particles is investigated using the 
inhomogeneous Balescu-Lenard equation expressed in terms of angle-action variables. The matrix method is implemented numerically 
in order to model the induced gravitational polarization. Special care is taken to account for the amplification of potential fluctuations 
of mutually resonant orbits and the unwinding of the induced swing amplified transients. Quantitative comparisons with A-body 
simulations yield consistent scalings with the number of particles and with the self-gravity of the disc: the fewer particles and the 
colder the disc, the faster the secular evolution. Secular evolution is driven by resonances, but does not depend on the initial phases 
of the disc. For a Mestel disc with Q~ 1.5, the polarization cloud around each star boosts up its secular effect by a factor of the order 
of a thousand or more, promoting accordingly the dynamical relevance of self-induced collisional secular evolution. The position and 
shape of the induced resonant ridge are found to be in very good agreement with the prediction of the Balescu-Lenard equation, which 
scales with the square of the susceptibility of the disc. 

In astrophysics, the inhomogeneous Balescu-Lenard equation may describe the secular diffusion of giant molecular clouds in galactic 
discs, the secular migration and segregation of planetesimals in proto-planetary discs, or even the long-term evolution of population 
of stars within the Galactic centre. It could be used as a valuable check of the accuracy of A-body integrators over secular timescales. 
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1. Introduction 


Galactic astronomy has striven to understand the dynamical evo¬ 
lution of discs over cosmic times. For these self-gravitating sys¬ 
tems, fluctuations of the potential induced by discrete encounters 
may be strongly amplified (Kalnajs 1972), while resonances tend 
to confine and localise their dissipation: such small stimuli can 
lead to long-term spontaneous evolution towards distinct galactic 
equilibria. The effect of this susceptibility on secular timescales 
will be addressed here in the context of an extended kinetic the¬ 
ory which takes explicitly into account such interactions. 

The kinetic theory of stellar systems was initiated by Jeans 
' 1929) and Chandrasekhar (1942) in the context of hot spherical 
stellar systems such as elliptical galaxies and globular clusters 
for which the gravitational susceptibility can safely be neglected. 
In contrast, self-gravitating galactic discs are cold dynamical sys¬ 
tems, for which rotation represents an important reservoir of free 
energy. More generally, in astrophysics, the secular diffusion 
of giant molecular clouds in galactic discs, the secular migra¬ 
tion and segregation of planetesimals in proto-planetary or de¬ 
bris discs, or even the long-term evolution of population of stars 
within the Galactic centre are all processes for which it is of 
interest to quantify the dynamical effect of gravitationally ampli¬ 
fied potential fluctuations induced by the finite number of stars 
involved. 

More than fifty five years ago, Balescu (1960) and Lenard 
(1960) developed a rigorous kinetic theory taking collective ef¬ 


fects into account, and obtained the corresponding kinetic equa¬ 
tion for plasmas, the Balescu-Lenard equation. More recently 
Heyvaerts (2010) and Chavanis (2012) have transposed the cor¬ 
responding non-linear kinetic equation to the angle-action vari¬ 
ables that are the appropriate variables to describe spatially in¬ 
homogeneous multi-periodic systems. The corresponding inho¬ 
mogeneous Balescu-Lenard equation accounts for self-driven or¬ 
bital secular diffusion of a self-gravitating system induced by 
the intrinsic shot noise due to its discreteness. Note that the for¬ 
mal transposition from position-velocity to angle-action implies 
that the secular interaction need not be local in space: they only 
need to correspond to gravitationally amplified long range corre¬ 
lations and resonances, which are indeed the driving mechanism 
for the secular evolution of isolated astrophysical discs via angu¬ 
lar momentum redistribution (Lynden-Bell & Kalnajs 1972). 


The Balescu-Lenard equation is valid at the order 1/A in 
an expansion of the dynamics in terms of this small parameter, 
where A» 1 is the number of stars. Therefore, it takes finite-A 
effects into account and describes the evolution of the system on 
a timescale of the order A/d, where to is the dynamical time. For 
self-gravitating systems, the collective effects are responsible for 
an anti-shielding which tends to increase the effective mass of 
the stars, hence reducing the relaxation time. When the system 
is cold, each particle is dressed by the very strong gravitational 
polarization it induces, hence the secular effects may occur on 
much shorter timescales than one would naively think, so that. 
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say iV e ff ~V/10 few . The purpose of this paper is to quantify this 
effect for stable but strongly susceptible galactic discs. 

The Balescu-Lenard formalism has seldomly been applied 
in its prime context, but only in various limits where it re¬ 
duces to simpler kinetic equations (Landau 1936; Vlasov 1938; 
Chandrasekhar 1942; Rosenbluth et al. 1957). Weinberg (1993) 
presents an interesting first implementation, though in a some¬ 
what over-simplified cartesian geometry. Yet, this formalism is 
quite unique in accounting for the non-linear evolution of discs 
and galaxies over secular timescales. V-body simulations, while 
potentially probing similar processes, should be scrutinized in 
such regime, as shadowing may, over many orbital times impact 
resonant interactions. V-body simulations have been shown 
to more or less reproduce growth rates of discs on dynamical 
timescales (see, e.g. Sellwood & Evans 2001, and references 
therein, together with Appendix C); qualifying them quantita¬ 
tively over secular timescales is now within reach of the Balescu- 
Lenard formalism. 

The companion paper, Fouvry et al. (2015a), hereafter pa¬ 
per I, presented a simple and tractable quadrature for the Balescu- 
Lenard drift and diffusion coefficients while assuming that the 
transient response of the disc was described by tightly wound 
spirals. Paper I applied the corresponding WKB approximation, 
while assuming that the disc was tepid and that the epicyclic 
approximation held. These simple expressions provided insight 
into the physical processes at work during the secular diffusion 
of self-gravitating discrete discs. When applied to the secular 
evolution of an isolated stationary self-gravitating Mestel disc, 
it identified the importance of the corotation resonance in the 
inner regions of the disc leading to a regime with both radial 
migration and heating, in qualitative agreement with numerical 
simulations. 

Yet, the tightly wound approximation is quantitatively ques¬ 
tionable when transient spirals unwind. Indeed paper I found 
a timescale discrepancy between the predicted secular evolu¬ 
tion timescale and the measured one, which might be driven 
by the incompleteness of the WKB basis. Such basis can only 
correctly represent tightly wound spirals. It also enforced local 
resonances, and did not allow for remote orbits to resonate, or 
wave packets to propagate between such non-local resonances. 
Yet, the seminal works from Goldreich & Lynden-Bell (1965); 
Julian & Toomre (1966); Toomre (1981) showed that any lead¬ 
ing spiral wave undergoes significant amplification during its un¬ 
winding to a trailing wave. Because it involves unwinding spirals 
this mechanism is not captured by the WKB formalism of paper 
I. 

In this paper, we will make no such approximations and will 
therefore compute numerically the corresponding diffusion and 
drift coefficients while relying on the matrix method (Kalnajs 
1976) to estimate the gravitational amplification of the secular re¬ 
sponse. It will allow us to assess the amplitude of the cross-talk 
between non-local resonances. We will then compare those pre¬ 
dictions to crafted sets of numerical experiments, allowing us to 
estimate ensemble averaged secular responses of a sizable num¬ 
ber of simulations as a function of the total number of particles V. 
Such ensemble average will allow us to make robust predictions 
for the V-scaling of the secular response and its dependence on 
halo to disc mass fraction, hence probing the secular importance 
of gravitational polarization. 

The paper is organized as follows. Section 2 briefly presents 
the content of the inhomogeneous Balescu-Lenard equation. Sec¬ 
tion 3 presents our implementation of the matrix method to com¬ 
pute the diffusion equation for an isolated self-gravitating ta¬ 
pered Mestel disc. Section computes numerically the exact 


drift and diffusion coefficient in action space for such a truncated 
Mestel disc, and compares the divergence of the corresponding 
flux density to the initial measured rate of change of the distri¬ 
bution function. Section presents our V-body simulations and 
compares scaling of the flux with the number of particles and 
the fraction of mass in the disc. Finally, section 6 wraps up. Ap¬ 
pendix A presents the relevant bi-orthogonal basis function. Ap¬ 
pendix C validates the response matrix method and the V-body 
integrator while matching growth rates and pattern speeds of un¬ 
stable Mestel discs. Appendix D investigates the roles of self¬ 
gravity and basis completeness. Appendix E describes the sam¬ 
pling strategy for the initial distribution. Appendix G presents 
briefly the available online codes. 


2. The inhomogeneous Balescu-Lenard equation 

We intend to describe the long-term evolution of a system made 
of V particles. We assume that the gravitational background il/o 
of the system is stationary and integrable, and associated with the 
Hamiltonian Ho- As a consequence, one can always remap the 
physical space-coordinates (x, v ) to the angle-action coordinates 
(0,J) (Goldstein 1950; Born 1960; Binney & Tremaine 2008). 
We define the intrinsic frequencies of motions along the action 
torus as 


Sl(J) = 0 = 


dH 0 

dj ' 


( 1 ) 


Within these new coordinates, one has that along the unper¬ 
turbed trajectories the angles 0 are 27r-periodic, evolving with 
the frequencies £2, whereas the actions J are conserved. We as¬ 
sume that the system is always in a virialised state, so that its 
distribution function (DF) can be written as a quasi-stationary 
DF of the form F-F{J,t), satisfying the normalization con¬ 
straint Jdxdv F = M tot , where M tot is the total mass of the sys¬ 
tem. On secular timescales, this isolated DF evolves under the 
effect of stellar encounters (finite-V effects). Such a collisional 
long-term evolution is descrided by the inhomogeneous Balescu- 
Lenard equation (Heyvaerts 2010; Chavanis 2012) which reads 


dF_ 

dt 


Tt(2Tt) d p 


dj i 


dj 2 


//ii ,mp ^ 


£ D (l?t r Qi-IM2-Q2) 
,mi (JuJ 2 ,m v £li)\ 
d 

dj 2 


, ( 2 ) 


where l/D miJ „,(Ji,J 2 , w) are the dressed susceptibility coeffi¬ 
cients, d is the dimension of the physical space, p = M tot /V is the 
mass of the individual particles, and where we used the short¬ 
ened notation = Since it is written as the divergence 

of a flux, this diffusion equation conserves the number of stars. 
One should also note the resonance condition encapsulated in 
the Dirac delta c>£,(mr£li-m 2 -n 2 ), with the integration over the 
dummy variable J 2 scanning for points where the resonance con¬ 
dition is satisfied. Note importantly right away that equation (2) 
scales like 1 /(ND 2 ) (since pccl/N), so that increasing V or in¬ 
creasing the heat content of the disc have the same effect. For a 
more detailed discussion on the content of the Balescu-Lenard 
equation, see paper I. 

In order to solve the non-local Poisson equation, we follow 
Kalnajs’ matrix method (Kalnajs 1976), so that we introduce a 
complete biorthonormal basis of potentials and densities i// (p Hx) 
and p (pl (x) such that 


At= 4 nGp (l ^ , 




dx [ijj (p \x)Y p (q \x) = -6 q p . 


(3) 
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The dressed susceptibility coefficients appearing in equation (2) 
are then given by 


1 

,//l2 


2 <Am,t/l) [I-M(W)]^ [^J 2 )T , 


(4) 


where I is the identity matrix and M is the response matrix de¬ 
fined as 

M W (W) = an) d Y f ij m dF/d J . (5) 

u co- mil 

m 

In the previous expression, we introduced as <p\n(J) the Fourier 
transform in angles of the basis elements if/ (p \x) using the con¬ 
vention that the Fourier transform of a function X(0, J ) is given 
by 

X mU)e im \ 

me f" r (6) 

X ^ = ^yJ^X(0,J)e— 0 . 


M from equation (5). Indeed, as noted in equation (3), the ma¬ 
trix relies on potential basis elements ([/ >pl which must be inte¬ 
grated over the whole action space with functions possessing 
a pole 1 /(ca-miA). This cumbersome and difficult evaluation 
has to be performed numerically, along with the matrix inver¬ 
sion needed to estimate the susceptibility coefficients from equa¬ 
tion (4). Finally, the third difficulty arises from the resonance 
condition fit which requires to determine how 

orbits may resonate one with another. In contrast, in paper I, we 
relied on the epicyclic approximation to build the angle-action 
mapping and on a WKB basis to treat gravity locally in order to 
solve these issues while obtaining tractable though approximate 
expressions. 

For a 2D axisymmetric potential, one can define explic¬ 
itly the actions of the system. Following Lynden-Bell & Kalnajs 
(1972); Tremaine & Weinberg (1984), the two natural actions of 
the system are given by a quadrature and an identity 

Ji=J r =-f dr j2(E-Mr))~L 2 /r\ 

n Jr p V (12) 

Ji — J<p = L, 


In order to ease the understanding of the Balescu-Lenard 
equation (2), one may rewrite it under the form of an anisotropic 
Fokker-Planck equation, by introducing the relevant drift and dif¬ 
fusion coefficients. Indeed, equation (2) can be put under the 
form 


dF_ 

dt 


y —■ 

4^ dJi 


, dF 

|A mi (/,)F(/i) + D mi (Ji)m v — 


(7) 


where r p and r a are respectively the pericentre and the apocen- 
tre of the trajectory, while E and L are the energy and angular 
momentum of the star. The first action /i encodes the amount 
of radial energy of the star, so that J r — 0 corresponds to circular 
orbits. The second action Jn is the angular momentum L of the 
star. One can then define the two intrinsic frequencies of motion 
Oi -k associated with the radial oscillations and TP = Tp, associ¬ 
ated with the azimuthal oscillations. Indeed, one has 


where A mi (J\) and D ltl] (J i) are respectively the drift and dif¬ 
fusion coefficients associated with a given resonance m. One 
should note that they both depend secularly on the distribution 
function, but this dependence was not exactly written out in or¬ 
der to shorten the notations. The drift coefficients are given by 


A mi (J) = -n(2n) d iJ 


?/ d 


6 D (/Hi-£Ji-/«2-^2) dF 

d/2—-—— - FT^ m2 7rr’ (8) 


and the diffusion coefficients are given by 


— = 2 f - 

Jr„ 


dr 


^2(E-Ur))-Jl/r 2 


03) 


while the azimuthal frequency Q 2 can then be determined via the 
relation 


O2 = h f'_ _dr_ 

Ql r 2 ^2(E-ip {) (r))-J 2 Jr 2 


(14) 


D mi (J l )=n(2n) a n) \dj 2 


?/ d 




;F(J 2 ) ■ (9) 


One can also introduce the total flux of diffusion ‘/ 7 tot as 

T m = Y J m [AmU)FU)+D m {J)m~ ) , 


( 10 ) 


so that the Balescu-Lenard equation from equations (2) and (7) 
takes the explicitly conservative form 


dF 

— = div CF tot ) . 
at 


(ID 


3. The Matrix diffusion equation 

When computing the Balescu-Lenard diffusion and drift coeffi¬ 
cients, three main difficulties have to be addressed. First, one 
must build the mapping (x, o)t-» (6, J), because the drift and 
diffusion coefficients are associated with a diffusion in action 
space. The second difficulty follows from the non-locality of 
Poisson’s equation and the estimation of the response matrix 


At this stage, one should note that various coordinates 
can be used to represent the 2D action space. Indeed, once 
the background potential ipo is known, one has the bijections 
(r p , r a ) <-> (E, L) <-2 (J r , J,p). As a consequence, any orbit can 
equivalently be represented by the set of the pericentre and apoc- 
entre (r p , r„) or by its actions (/1 , Ji). However, determining the 
actions associated with one set ( r p , r a ) only requires the compu¬ 
tation of a ID integral as in equation (12), whereas determing 
the pericentre and apocentre associated with a set of actions 
requires the inversion of the same non-trivial relation. 
Because the peri/apocentres are the two roots of the equation 
2(E-if/o(r))-L 2 /r 2 =0, one also immediately obtains that for a 
given value of r p and r a , the energy E and the angular momen¬ 
tum L of the orbit are immediately given by 

r 2 ifi a — r 2 (f/ p \2(iJf a -<Jj p ) 

F r 2 — r- ’ ^ V r- 2 -r- 2 ’ 

'a 'p y ' p 'a 

where we used the shortening notations i//p/a-^o( p p/a)- There¬ 
fore, in the upcoming calculations, we will use ( r p , r a ) as the 
representative variables of the 2D action space. 
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3.1. The basis elements 


3.2. Computation of the response matrix 


The expressions (4) and (5) of the susceptibility coefficients and 
the response matrix require the introduction of 2D potential- 
density basis elements. The 2D potential basis elements that 
we will consider will depend on two indices spanning the two 
degrees of freedom so that one has 

^\R, b) = ti(R, <P) = e lt<p <U ( n {R, cf> ), (16) 

where T( f H is a real radial function and (R. <p) are the usual polar 
coordinates. The associated surface densities elements will be of 
the form 


z {p) (R, <f>) = K(R, 0) - e >e * &n(R, <f >), (17) 

where fD e n is a real radial function. The basis elements therefore 
depend on two indices £>0 and n> 0. In all the numerical calcu¬ 
lations, we used the radial functions from Kalnajs (1976), which 
are recalled in Appendix A. 

The next step is then to determine the Fourier transform with 
respect to the angles 6 of the basis elements. Indeed, to compute 
the response matrix M from equation (5), one has to compute 
where the resonance vector is given by m - (ni\, m 2 ). Fol¬ 
lowing equation (6), it is given by 

f iSidffo «A W (R, 4>) e~ im ' e ' e~ imA . (18) 

From Lynden-Bell & Kalnajs (1972), the angles 8\ and 8i asso¬ 
ciated with the actions from equation (12) are given by 


fli = n. 

[dr 1 

,c > ^2{E-Mr))-J*/r 

+ 15 
■©- 
11 

OS 

S'-^ 

Th-h/r 2 

- ^ 2(E-Mr))~Jj/r 2 


(19) 


where Cj is a contour starting from the pericentre r p and going 
up to the current position r = r(8\) along the radial oscillation. 
Following the notations from Tremaine & Weinberg (1984), one 
can straightforwardly show that equation (18) takes the form 

^(J) = 8Zw;; m2 n P (J), (20) 

where f W^ m2nP U) is given by 

ir (r) cos[m , 0, [r]+m 2 (0>-0)[r ]]. (21) 

xJr P 

In equation (21), the boundaries of the integral are given by the 
pericentre r p and apocentre r„ associated with the action J. Such 
an expression underlines the reason why ( r p , r a ) appear naturally 
as good coordinates to describe the 2D action space. One can 
note that equation (21) involves an integral over r thanks to the 
change of variables 8\ —> r, which satisfies 


d0i _ Qi 

dr ^2(E-Mr))~Jl/r 2 


( 22 ) 


In equation (21), 8\[r] and (02 — 0)[^] only depend on r via 
the mappings from equation (19). Provided that r U pP is a real 
function, the coefficients < TVjp m2nP are always real. Because 
these coefficients involve two indicated integrals, they are nu¬ 
merically expensive to compute. However by parity, they obey 
= 'W'fj'mjnr ’ which allows a significant reduction of the 
number of coefficients to compute. 


We now have all the elements required to compute the response 
matrix from equation (5). In its definition, one should note the 
presence of an integral over the mute variable J , which, as dis¬ 
cussed previously, will be performed in the 2D ( r p , r fl )-space. 
The first step is to go from J = (J\ , J 2 ) to (E, L). The Jacobian of 
this transformation is given by 


d(E,L) 

d(Ji,J 2 ) 


dE 

dE 



dJ\ 

~dT> 


Qj O 2 

dL 

dL 


0 1 

dJ\ 

dh 




(23) 


so that one immediately has d/id /2 =d£'dL /El\. Given the ex¬ 
pression (20) of the 2D Fourier transformed basis elements, the 
response matrix may be written under the form 


M p Jaj) = (2nY6] 


2 ^Y [dEdL^- 


m [ * 


Qj 


(m u {P)-dFo/dJ 


xW 


m 1 

tP£PnP 




(24) 


where the sum on m 2 has been dropped. Moreover, we dropped 
the conjugate over 'W m [p[PnP since they are always real. We may 
now perform the change of variables (E, L)—*{r p ,r a ), so as to 
rewrite equation (24) under the form 


M p Jto) = 




JPnVrfli 


dr p dr. 


9m " OpDa) 
K mtP (r p ,r a ) ’ 


(25) 


where the functions g l lp ' Pn \r p , r a ) and h“ (p (r p ,r a ) are respec¬ 
tively given by 


gZ? n \r p ,r a )=(2n) 2 


d(E, L ) 

1 

d(r p , r a ) 





dF 0 

dj 




and 


(26) 


K lieP {r p ,r a ) = io-{m l J p )-£l. (27) 

It is important to note here that the response matrix is diago¬ 
nal with respect to the ( p and ( q indices so that each t may 
be treated independently. The definition of the function g from 
equation (26) involves the Jacobian d(E, L)/d(r p , r a ) of the trans¬ 
formation (E, L) —» (r p , r a ) which can be immediately computed 
from the expressions (15) of E-E{r p , r a ) and L-L(r p , r a ). More¬ 
over, in some situations, the DF F = F(J), may also rather be de¬ 
fined as F -F(E,L). It is straightforward to show that one has 


dF „ dF 

m sJ‘ m ' a, aE 



' dF 

dF 


+ m 2 

L 

^ 2 dE 

l + ~dE 

E 


(28) 


3.3. Sub-region integration 

The next step of the calculation is then to perform the remain¬ 
ing integration over ( r p , r„ ) from equation (25). Because of the 
presence of the resonant pole 1 //?" (p , such a numerical integra¬ 
tion has to be performed carefully. We cut out the integration do¬ 
main (r p , r„) in various subregions indexed by i. The /^-region 
will be centred around the position ( F , r‘ a ) and will correspond 
to the square domain such that r p e\r' p -Ar/2 ; r' 7 + Ar/2] and 
r a e[r' a -Ar/2\ r^+Ar/2], where A r corresponds to the size of 
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the subregions. The smaller A r, the better will be the approxi¬ 
mated estimations of the response matrix. Within the I th — region, 
one can write first-order Taylor expansions of the functions g and 
h from equations (26) and (27) around the centre (r' p , r‘ a ) of the 
region such that 

fg(r‘ p + Ar p , r' a +Ar a ) ^ a l g +b l g Ar p +d g Ar a , 

\li(r' p +Ar p , r' a +Ar a ) a\+b\Ar p +c' h Ar a , 

where for convenience we shortened the index dependences from 
equations (26) and (27). The coefficients a‘ g , b g and c' g (similarly 
for h) are given by 


a' =g(r;,r'); b‘ g = jjj- 


(d,r‘a) 


C 9 = 


dg 

dr a 


(A,A) 


(30) 


3.4. Critical resonant line 

The resonance condition encapsulated in the Dirac delta 
6 d(uii £ii -m 2 -Hi) generates an additional difficulty in the 
calculation of the Balescu-Lenard drift and diffusion coefficients 
from equations (8) and (9). Recall the definition of the composi¬ 
tion of a Dirac delta and a function (Hormander 2003), which in 
a af-dimensional setup takes the form 

f dxf(x)S D (g(x)) = f dcr(x) _ /U * , (34) 

Jg-'i 0) \Vg(X)\ 

where g _ 1 ( 0 ) = {jt | g(x)- 0 ) is the hyper-surface of dimension 
(generically) id- 1) defined by the constraint g(x)- 0 , and d cr(x) 
is the surface measure on 0 _ 1 (O). We have also defined |Vgf| as 
the euclidean norm of the gradient of g , so that one has 


where it is important to note that these coefficients are only func¬ 
tions of the central coordinates ( r' p , r' u ) and will be treated as con¬ 
stants on each sub-region. In the numerical implementation, the 
coefficients involving partial derivatives will be estimated by fi¬ 
nite differences, so that one will have for instance 


bg(r‘ p ,r l a ) 


gidp+Ar.dJ-gidp-Ar,^) 

2Ar 


(31) 


which allows to minimize the number of evaluations of g re¬ 
quired. The approximated integration on each sub-region can 
then be performed and takes the form 


lit 


dr p dr a 


g(r p , r a ) 
h(r p , r a ) 


d.r„d.r„ 


a'g + b'gXp+c l g x a 

a l h +b l h x p +c l h x a +ir] 


= N(a', b‘ c' a‘ h , b l h , c‘ h , rj, Ar), 


(32) 


where N is an analytical function which only depends on the coef¬ 
ficients obtained in the limited developments from equation (29). 
In order to have a well-defined integral, we added an imaginary 
part 77 > 0 to the temporal frequency w, so that a>-a>o + ig. When 
looking for unstable modes in a disc, this imaginary part 77 corre¬ 
sponds to the growth rate of the mode, ft is also crucial to note 
here that one always has a' g , b' g , c^eR and similarly a ' h , b' h , c^eR. 
The effective computation of the function N is presented in Ap¬ 
pendix B. Thanks to equation (32), the expression (25) becomes 


M pq (aS) = Yj s H’ b r c r a ‘h - b h> c 'h - 7, Ar). (33) 

mj i 

In the previous expression, in order to effectively compute nu¬ 
merically the sum on m\, we introduce a bound m” ax , so that the 
sum is only reduced to |mi|<m™ ax . Because of the requirement 
to truncate the action space in various subregions as in equa¬ 
tion (33), the computation of the response matrix still remains 
a daunting task, to ensure appropriate numerical convergence. 
In Appendix C, we detail the validation of our implementation 
of the response matrix calculation, by recovering known unsta¬ 
ble modes of truncated Mestel discs (Zang 1976; Evans & Read 
1998b; Sellwood & Evans 2001). Once the response matrix M 
is known, the determination of the dressed susceptibility coeffi¬ 
cients 1/|1D | 2 from equation ( !) involves a straightforward sum¬ 
mation . 


1 One could if needed regularize the inversion of I-M, to avoid Gibbs 
rigging, since our basis is significantly truncated; this has proven not 
necessary here. 


|V<7(jc)| = 


dg 

2 

dg 

i dxi 

+ ...+ 

dx d 


(35) 


Here we have assumed that the resonance condition associated 
with the function = is non-degenerate , so 

that Vjteg _ 1 (0), |Vg(jc)| >0, which also ensures that the dimen¬ 
sion of g 1 (0) is (d- 1). One should note that this degeneracy 
condition is not satisfied by the Keplerian or harmonic potentials. 
Because we are considering an infinitely thin disc, the dimension 
of the physical space is given by d- 2 , so that the set g 1 ( 0 ) will 
take the form of a curve y, that we will call the critical resonant 
curve. Generically, it will take the form of an application of the 
type 


y : M6[0;l]Hy(tt) = (7 1 (#)j 2 («)). 


(36) 


One can then immediately rewrite the r.h.s of equation (34) under 
the form 


i 


dcr(.t) 


fix) 

|v<Kjc)| 


■X 


= I dw 


fiyiuj) 

|v<Kr(H))l 


lr'(«)l, 


(37) 


where we have naturally defined \y'(u)\ as 

\y\u)\ = 


dyi 


du 


dyi 


du 


(38) 


Therefore, as soon as the critical resonant curve y has been iden¬ 
tified, the integration from equation (34) can be computed. 

As noted in equation (15), using the pericentres and apocen- 
tres ( r p , r a ). given the Jacobian from equation (23) and proceed¬ 
ing in the same way as in equation (25) for the response matrix, 
one may rewrite the drift and diffusion coefficients from equa¬ 
tions (8) and (9) under the form 

A mi (J 1 ) =Y fdr p dr a 6 D (m v £l l -m 2 -£l 2 ) Gf ni „, 2 (r p , r a ), (39) 

mj ^ 

and 

D mi (Ji) =Y f'drpdr a S D (mr£li-m 2 -Sl 2 )G° 1 mi (r p ,r a ). (40) 

mj ^ 


where the functions G,^ w and G|,’ mi are respectively defined 
as 


ci 


Xr p ,r a ) = - 


1 

Q[ 


d(E.L) 4n 3 gtti2-dF/dJ2 

d(r p , r a ) |D m , ,m 2 (/i, Ji, m -^i)l 2 


(41) 
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and 


r D 

''- 7 /w l .nil 


(r p , r a ) 


1 

Q) 


d(E,L) 4ApF{Ji) 

d(r p ,r a ) \D mi , m2 (J u J 2 ,m v Sl l )\ 2 ' 


(42) 


For a given value of J \ , m i and mi , and introducing a>i = m\ ■ il i, 
we define the critical curve y mi (<x> i) as 


y,n 2 («i) = |(r p ,r fl )|in 2 -il(r / „r fl ) = (yi|. (43) 

The expressions (39) and (40) of the drift and diffusion coeffi¬ 
cients immediately become 


^4/wi (J 1) 
D mi (J \) 


-Zf 


/~<A 

^ni\ ,ni2 


= 2 / 


dcr ■ 

y^(«i) |V(iw 2 -£i2)| ' 

r D 

, ^//li ,/Hi 

dcr-— 

Ym 2 (0>l) |V(//l2'i^2)l 


(44) 


4.1. Initial setup 

We consider the same disc as considered in Sellwood (2012). It 
is an infinitely thin Mestel disc for which the circular speed v<p 
is a constant Vo independent of the radius. The stationary back¬ 
ground potential lAyi and its associated surface density I\i are 
given by 


if/u(R) = V,] log 


R 

Rmax 


Zm(R) - 


Vq 

2kGR ' 


(46) 


where R max is a scale parameter of the disc. Following Toomre 
(1977); Binney & Tremaine (2008), a self-consistent DF for this 
system is given by 

F m (E, J 0 ) = C M Jl exp[-£/cr;], (47) 

where the exponent q is given by 


where the resonant contribution |V(wi 2 '^ 2 )[ is defined as 


|V(/w 2 -fi 2 )l 


d il 2 

2 

d£l 2 

W2-—- 

dr p 

+ 

m 2 ■ —— 

dr a J 


(45) 


The derivatives of the intrinsic frequencies with respect to r p 
and r a appearing in equation (45) will be computed as in equa¬ 
tion (3]) using finite differences. Once the critical lines of reso¬ 
nances have been determined, the computation of the drift and 
diffusion coefficients from equation (44) is straightforward, so 
that the full secular diffusion flux T XM from equation (10) may 
be determined. 




1, 


(48) 


with cr r being the constant radial velocities spread within the 
disc. In equation ( 17), Cm is a normalization constant given by 


Cm - 


V 2 

v o 


2 l +? /2 7r 3/2 Gcr ?+ 2 r [l + |] jR ? + a l 


(49) 


In order to deal with the central singularity of the Mestel disc 
along with its infinite extent, we introduce two tapering functions 


4. Predicting Balescu-Lenard flux divergences 

We may now illustrate how the previous computations of the 
response matrix and the Balescu-Lenard drift and diffusion 
coefficients can be used to recover some results obtained in 
well-crafted numerical simulations of galactic discs. Indeed, 
Sellwood (2012) (hereafter S12), studied the long-term evolution 
of an isolated stable truncated Mestel disc (Mestel 1963). After 
letting the disc evolve for hundreds of dynamical times, S12 ob¬ 
served a secular diffusion of the disc DF in action space, through 
the spontaneous generation of transient spiral waves. The most 
striking result of this evolution is given in figure 7 of SI 2, which 
exhibits the late time formation of a resonant ridge in action 
space along a specific resonant direction. Such diffusion features 
observed in the late evolution of an isolated stable and discrete 
system are thought to be signatures of a secular evolution in¬ 
duced by fi n ite—/V effects, as described by the Balescu-Lenard 
formalism. Because the system is made of a finite number N of 
pointwise particles, it undergoes (long range) resonant encoun¬ 
ters leading to an irreversible secular evolution. In order to in¬ 
vestigate such a collisional evolution, paper I applied the WKB 
limit of the Balescu-Lenard formalism to S12 simulation. While 
most of the secular diffusion was qualitatively recovered, there 
remained a significant timescale discrepancy, since the typical 
timescale diffusion predicted by this approach was typically 10 3 
times too slow compared to the observations made in S12. The 
use of a non-local basis such as equation ( ) and the numerical 

computation of the response matrix from equation (24) allows to 
incorporate in the present paper these previously ignored contri¬ 
butions from the WKB approach. In the upcoming sections, we 
therefore present briefly the disc considered by S12 and our de¬ 
termination of the secular diffusion flux predicted by the Balescu- 
Lenard formalism. 


J Vl 

T’inner(^) = * 

7outcr(7 (f>) — 

where the indices v t and p t control the sharpness of the two ta¬ 
pers, while the radii R x and Rq are two scale parameters. These 
tapers 7) nncl and T outer respectively represent the bulge and the 
outer truncation of the disc. In addition to these taperings, we 
also suppose that only a fraction £ of the stellar disc is self- 
gravitating, with 0<£ < 1, while the rest of the gravitational po¬ 
tential is provided by the static halo. As a consequence, the active 
distribution function F star is given by 

F’star(£’, Jcj,) = £ Fm(E, J,/,) -dinner tp) ^outer (A)- (51) 

We place ourselves in the same units system as in S12, so that 
we have To =G = R, = \. The other numerical factors are given 
by #=11.4, v t =4, p t -5, £ = 0.5, /?o = H-5 and /? max = 20. The 
contours of the tapered DF F sXal - are illustrated in figure . At 
this stage, it is important to note that S12 restricted the pertur¬ 
bations forces to the harmonic sector m<p = 2, so that we may 
consider the same restriction on the considered azimuthal num¬ 
ber m ,p. As a consequence, in the double resonance sum on m\ 
and m 2 present in the Balescu-Lenard flux from equation ( 2 ), 
we will assume that m\ and mi belong to the restricted set 
{m r ,m ^}G{(-1,2),(0,2),(1,2)}, where — 1,2) corre¬ 

sponds to the Inner Lindblad resonance (ILR), (m r ,m^)-( 0,2) 
to the Corotation resonance (COR), and ( m r , m,f) = (1,2) to the 
Outer Lindblad Resonance (OLR). All the upcoming calcula¬ 
tions have also been performed while taking into account the con¬ 
tributions from the resonances associated with m, = ±2, which 
were checked to be largely subdominant. 


AVq^+j; 
Jd, 


(50) 
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2 4 6 8 10 

dip 

Fig. 1: Contours of the initial distribution function F sta r from equa¬ 
tion ( 51 ), in action space The contours are spaced linearly be¬ 

tween 95% and 5% of the distribution function maximum. 


4.2. Initial drift and diffusion 

As detailed in equation (33), the computation of the response 
matrix requires to consider a grid in the ( r p , r„)- space. We con¬ 
sidered a grid such that r™ n =0.08, r™ ax = 4.92 and Ar = 0.05. 
The sum on ni\ appearing in equation (33) was reduced 
to |mi|<m” ax = 7. The basis considered was Kalnajs 2D ba¬ 
sis (Kalnajs 1976) with the parameters k Kd = l and R^-5. One 
should note that despite having a disc which extends up to 
R max = 20, one can still consider a basis truncated at such a small 
R k : i , so as to be able to efficiently capture the diffusion proper¬ 
ties of the system in its inner regions, from where the secular 
diffusion is known to start. The radial basis elements were re¬ 
stricted to 0<«<8. When evaluating the response matrix, as in 
equation (32), one has to add a small imaginary part r/ to the 
frequency so as to regularize the resonant denominator. Through¬ 
out the calculations presented below, we considered r/= 1CT 4 and 
checked that this choice had no impact on our results. 

Since the total potential i^m is known via equation ( 16), 
the mapping to the angle-action coordinates is completely de¬ 
termined. The two intrinsic frequencies of the system can then 
be computed on the (r p ,r a )-grid via equations ( 3) and ( ). 

Once these frequencies are known, the critical resonant lines in¬ 
troduced in equation (43) can be determined and are illustrated 
in figure 2. It is along these lines that one will have to perform 
the integration present in the definitions of the drift and diffusion 
coefficients from equation (44). 

Thanks to this expression, one can then compute the secular 
diffusion flux Tun defined in equation ( ). Because the mass 

of the particles is given by /j-M tot /N, it is natural to consider 
the quantity NT tot which is independent of N. The vector field 
-NTmi - -(N'Fj ir NTj r ), which represents the direction of dif¬ 
fusion of individual particles, is illustrated in figure 3. One can 
already note in figure 3 that the diffusion vector field is along 
a narrow resonant direction. Along this ridge, one typically has 
Tjf——2Tj r , so that the diffusion appears as aligned with the di¬ 
rection of the ILR resonance given by /«ilr = (2, -1). If one con¬ 
siders only the curl free part of this vector field, a sink and a 
source can be easily identified within that flow. 

Once the diffusion flux NTtot has been determined, one can 
compute the divergence of this flux, so as to determine the re- 



Fig. 2: Illustration of 4 different critical resonant lines in the 
(r p , r„)-space. As defined in equation (43), a critical line is character¬ 
ized by the resonant vectors m t , m 2 and a location J\ <-> (r' p , R a ) in ac¬ 
tion space. Each of the 4 plotted critical lines are associated with the 
same location (r p , rj), represented by the black dot. The critical lines 
correspond to various choices of resonant vectors mi and m 2 among the 
three inner outer and corotation Lindblad resonances, respectively ILR, 
OLR and COR. One should also note that for m\ = m 2 , the critical lines 
go through the point (r p , r*). 



Fig. 3: Map of -NTtot, where the flux has been computed with 
mi, m 2 e{mi L R, m c or, /»olr1- Following equation (1 ), —NT tot corre¬ 
sponds to the direction of diffusion of individual particles in action- 
space. 


gions for which the DF is expected to change during the secu¬ 
lar diffusion. Figure illustrates the contours of A r div(y r tot ). In 
figure 4, we obtained that the Balescu-Lenard formalism indu¬ 
bitably predicts the formation of a narrow resonant ridge aligned 
with the ILR-direction, as was observed in S12 simulation. One 
also recovers that the stars which will populate the resonant ridge 
originate from the basis of the ridge and diffuse along the direc¬ 
tion associated with the ILR resonance. It is most likely that the 
slight shift in the position of the ridge is due to the fact that 
the Balescu-Lenard prediction was carried at t = 0 + , while S12’s 
measurements are at t - 1400, so that we do not expect a perfect 
match. Other sources of discrepancies might be the use of a soft- 
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Fig. 4: Left panel : Map of AdivCFtot), where the total flux has been computed with m { , m 2 E{m lLR , m C oR, »Iolr}- Red contours, for which 
iVdiv(7 r tot)<0 are associated with regions from which the orbits will be depleted, whereas blue contours, for which AdivCFtot)>0 correspond 
to regions where the value of the DF will be increased during the secular diffusion. The contours are spaced linearly between the minimum and 
the maximum of Miiv(Ftot)- The maximum value for the positive blue contours corresponds to MIivOFtot)-350, while the mininum value for 
the negative red contours is associated with Adiv(;F t ot)- _ 250. Right panel : From Sellwood (2012) - figure 7, contours of the change in the DF 
between the time fsi 2 = 1400 and fsi 2 = 0, for a run with 50 M particles. Similarly to the left panel, red contours correspond to negative differences, 
i.e. regions emptied from their orbits, while blue contours correspond to positive differences, i.e. regions where the DF has increased during the 
diffusion. Both of these contours are aligned with the ILR direction of /n 1LR = (2, -1) in the (7^, 7 r )-plane, corresponding to the cyan line. 


ening length in numerical simulations, which modifies the two- 
body interaction potential, or the difference between an ensem¬ 
ble average (as predicted by the secular formalism) and one spe¬ 
cific realization - our own simulations suggest that there is some 
variation in the position of the ridge between one run and an¬ 
other. Because we explicitly determined the value of NdivCF tot ), 
we may now study the typical timescale of collisional relaxation 
predicted by this Balescu-Lenard estimation as detailed in sec¬ 
tion -.3. One may also investigate the respective roles of the self- 
gravitating amplification and the limitation to the tightly-wound 
basis elements as presented in Appendices and D.2. 

4.3. Timescale of diffusion 

The most significant disagreement found in paper I, while apply¬ 
ing the WKB approximation of the Balescu-Lenard equation to 
S12’s simulation was a timescale discrepancy between the time 
required to observe the resonant ridge in S12 simulation and the 
collisional timescale for which the finite-/V effects come into 
play. As already noted in paper I, since the Balescu-Lenard equa¬ 
tion (2) only depends on N through the mass of the individual 
particles p-M tot /N, we may rewrite it under the form 

(52, 

where Cbi.[F] = AdivFFtot) is the A'-indcpendent Balescu- 
Lenard collisional operator, i.e. the r.h.s of equation (2) mul¬ 
tiplied by N — M tot /p. As expected, the larger the number of 
particles, the slower the secular evolution. This also illustrates 
the fact that the Balescu-Lenard equation comes from a kinetic 
Taylor expansion in the small parameter s-l/N^l. Intro¬ 
ducing the rescaled time r-t/N , so that equation ( 2) reads 

C\ jp 

= Cbl[F] , (53) 

letting us express the Balescu-Lenard equation without any ex¬ 
plicit appearance of N. In paper I, we estimated the time Arsn 


required to observe the ridge as Arsn — 3 x 1CL 5 . When perform¬ 
ing the same measurement thanks to the contours of the diffu¬ 
sion flux div(7 r tot ) computed within the WKB approximation, 
we obtained Atwkb — 3x KL 2 , so that paper I obtained the ra¬ 
tio Atsi 2 /Atwkb — 1CL 3 . This discrepancy was due to the lim¬ 
itation to tightly wound spirals. Because the estimation of the 
secular diffusion flux Tioi presented in figure 4 was made us¬ 
ing the matrix method (Kalnajs 1976) with a full basis, it cap¬ 
tures the additional swing-amplification. Indeed, given the map 
of Adiv'Ftot obtained in figure 1, one may estimate the typical 
time Atbl required for such a flux to lead to the diffusion features 
observed in S12. The contours presented in figure 7 of S12 are 
separated by an increment equal to 0.1 xfj“, where F” ax -0.12 
is the maximum of the normalized DF (via equation (51)). In or¬ 
der to observe the resonant ridge, the value of the DF should 
typically change by an amount of the order AFo^0.1 xF™. 
From figure 4, one can note that the maximum of the diver¬ 
gence of the diffusion flux is given by |Adiv(!Ftot)lmax-350. 
Thanks to equation ( ), one can immediately write the rela¬ 

tion AFo- ATBiJATliv(!Ftot)|max, where Atbl is the time during 
which the Balescu-Lenard equation has to be evolved in order 
to develop a ridge. With the previous numerical values, one ob¬ 
tains Atbl — 3 X 10 \ Comparing the numerically measured time 
Atsi 2 and the time Atbl predicted from the Balescu-Lenard 
equation, one obtains 


Atsi2 

At B l 


(54) 


As expected, the projection of the response over an unbiased ba¬ 
sis leads to over a hundredfold increase of the susceptibility of 
the disc and therefore to a very significant acceleration of sec¬ 
ular diffusion. Thanks to this mechanism, we now find a very 
good agreement between the diffusion timescales observed in nu¬ 
merical simulations and the predictions from the Balescu-Lenard 
formalism. This quantitative match is rewarding, both from the 
point of view of the accuracy of the integrator (symplecticity. 
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timestep size, softening...), and from the relevance of the succes¬ 
sive approximations underpinning the Balescu-Lenard formal¬ 
ism (timescale decoupling, truncation of the BBGKY hierarchy, 
neglect of the close encounter term...). 

In Appendix D, we show that when considering either fig¬ 
ure D. 1, for which the self-gravity of the system has been turned 
off, or figure D.3 for which the loosely wound basis elements 
were not taken into account, one does not recover a narrow reso¬ 
nant ridge appearing on timescales compatible with S12 simula¬ 
tions. Therefore, the main source of secular collisional diffusion 
oberved in S12 and recovered in figure has to be the strong 
self-gravitating amplification of loosely wound perturbations, i.e. 
a sequence of uncorrelated swing-amplified spirals sourced by 
finite-/V effects is indeed the main driver of secular diffusion. 
The WKB formalism from paper I identified correctly the fam¬ 
ily of orbits involved, but fell short in predicting how narrow the 
resonant ridge is and how strongly amplified the response is. 

5. Comparison to N -body simulations 

In paper I, we relied on simulations presented in Sellwood (2012) 
to compare the divergence of the diffusion flux to the WKB pre¬ 
diction. In order to probe the expected scalings with the number 
of particles or with the active fraction of the disc, we now resort 
to our own /V-body simulations. 

5.1. N-body integration 

The initial sampling of particles is critical when investigating 
the origin of secular evolution, as one must ensure that the disc 
is initially in a state of equilibrium. The sampling strategy we 
implemented is described in some detail in Appendix . 

Once sampled, we evolve the initial conditions using a 
straightforward particle-mesh A'-body code with a single- 
timestep leapfrog integrator (e.g., Binney & Tremaine 2008, 
§3.4.1). We follow S12 and split the potential in which the 
particles move into two parts: (i) an axisymmetric contribution 
ijj m from the unperturbed Mestel disc, as in equation ( 6) and 
(ii) a non-axisymmetric contribution \p\(R,cp) that develops as 
perturbations grow in the disc. This splitting avoids difficulties 
in the treatment of the rigid component of the potential that 
is not included in the DF, due to the tapering functions and 
active fraction introduced in equation (51). We calculate t/q 
using cloud-in-cell interpolation (e.g., Binney & Tremaine 2008, 
§2.9.3) of the particles’ masses onto an A mes h X A mes h mesh of 
square cells spaced Ax apart, then filtering the resulting density 
held p(x,y) to isolate the disc response (see below), before ap¬ 
plying the usual Fourier-space doubling up procedure to obtain 
the potential t/q at the cell vertices. The contribution of t/q to 
each particle’s acceleration is then obtained using the same 
cloud-in-cell interpolation scheme. 

When computing the density mesh, we added a filtering 
scheme, to include only the m^-2 disc response, similarly to 
what was considered in S12. We isolate this m$ = 2 mode by cal¬ 
culating 

Pi(r) = ^-Jd0p(rcos(<^),rsin(0))e~ 2 ^, (55) 

immediately after the cloud-in-cell assignment of mass to the 
mesh at each timestep, then imposing the new mesh mass distri¬ 
bution 

p(xk, yk ) = Piin) e 2 ' 0t , (56) 


with (r/., (pi,) chosen according to (**, yk) = (n cos(0*), r* sin(0£)). 
To obtain pi(r) we use brute-force computation of equation (55) 
on a serie of A' Ilng radial rings with spacing Ar*szAx, using the 
trapezium rule with N& =720 points in </> for the angular integrals. 
These models are designed to reproduce as closely as possible 
the essential details of S12’s simulations. There are a couple of 
deliberate technical differences: S12 uses a polar mesh to obtain 
i/q, whereas we use a cartesian mesh with am^=2 prefiltering of 
the density held; S12 has a block timestep scheme instead of our 
simpler single-timestep one. 

For the results presented here we used a timestep 
At = 10 3 /?i/Vo on a mesh that extends to ±f? max = 20with 
A^mesh = 120 cells, so that Ax=R J3. The filtering of the potential 
perturbations to the harmonic sector m^- 2, was performed with 
Wring = 1000 radial rings, so that Ar-2 R;/100, and A^ = 720 
points in the azimuthal direction. Finally, the computation of the 
potential from the density via Fourier transform, was performed 
with a softening length s-R{/ 6, which is comparable to the 
value used in Sellwood (2012), which considered a Plummer 
softening with s-R{/ 8. The results are not significantly changed 
when we halve the timestep or the mesh size. In Appendix C, 
we detail the validation of our A-body code, by recovering 
known unstable modes of truncated Mestel discs (Zang 1976; 
Evans & Read 1998b; Sellwood & Evans 2001). 

5.2. Seating with N 

In order to rid our measurements of individual fluctuations, we 
run multiple simulations for the same number of particles and 
perform an ensemble average of different evolution realizations 
for the same number of particles. It allows us to estimate only 
the mean evolution , which is effectively what is described by the 
Balescu-Lenard formalism. 

In order to study the scaling with A of these numerical sim¬ 
ulations, one has to extract from the simulations a quantity on 
which to test this scaling and compare it with the predictions 
from the Balescu-Lenard formalism. The statistical nature of the 
initial sampling presents an additional difficulty. Indeed, because 
one only samples A stars as described in Appendix E, the initial 
effective DF fluctuates around the smooth background DF from 
equation (51) as a Poisson shot noise. These statistical fluctu¬ 
ations originate from the initial sampling and are not as such 
specific to the physical process captured by the Balescu-Lenard 
formalism, so that one should carefully disentangle these two 
contributions. Hence we introduce the function h(t,N ) defined 
as 

h(t,N) = (hi(t,N)) , (57) 

where the operator (■) corresponds to the ensemble average, ap¬ 
proximated here with the arithmetic average over the p = 32 dif¬ 
ferent realizations of simulations for the same number of parti¬ 
cles, indexed by i: (■)- 1 / pX, (■). In equation (57) the function 
hi(t, N) is a lag function which read 

hi(t,N) =Jd J [E,(f,y,A)-<E(f = 0,/,A))] 2 , (58) 

where we defined as F,(t,J,N) the normalized DF of the I th re¬ 
alization for a number A of particles. Such a quantity intends to 
quantify the distance between the initial mean DF (F(t- 0)) and 
the evolved DF /q. We are interested in the early time behavior 
of the lag function h from equation (57), so that we may perform 
its Taylor expansion 

~ t 2 

h(t , N) =* ho(N) + h\(N) t + h 2 (N) - , (59) 
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where it is important to note that the coefficients hg , h\ and h 2 
depend only on N and are given by 


h 0 (N) = h(t=0,N) ■ hi(N) = ^ 

at 


d 2 h 

; h 2 (N) = — 
t =0 ot- 


7=0 


(60) 


Let us now estimate each of these coefficients in turn. Thanks to 
equation (58), one can compute h q(N) which reads 



{[F-(Fg)] 2 ) , 


(61) 


where we used the shortened notations ( Fq ) = (F(t = 0, N)) and 
F-F(t-0,J,N). We note that this coefficient only depends on 
the properties of the initial sampling, and not on its dynam¬ 
ics. Because discrete sampling obeys Poisson statistics, one can 
write 

h 0 (N) = ^ , (62) 

where ag is a constant independent of N. One may then compute 
/zi(lV), which takes the form 

hm = 2Jdj ([F-(Fo)] F') , (63) 


where we used the shortened notation F' = [dF/dt](t= 0). One 
should note that the terms appearing in equation (63) have two 
different physical contents. Indeed, the term [T'-(T’o)] involves 
initial sampling, whereas F' is driven by the dynamics of the sys¬ 
tem. If we assume that the sampling and the system’s dynamics 
are uncorrelated, one writes 


{[F-(F 0 )]F') = (F-(F 0 ))(F')=0. 


(64) 


As a consequence, one immediately obtains from equation (63) 
that h\(N) = 0. One can finally compute the coefficient JioiN) 
which reads 

h 2 (N) = 2jdJ ([F'] 2 + [F-(F 0 )1 F") , (65) 

where we used the shortened notation F” = [d 2 F/dt 2 \(t = 0). Us¬ 
ing the same argument as in equation (6‘i), we may get rid of the 
second term in the l.h.s of equation (65). If we now also assume 
that the variance of [ F'] 2 is small compared to its expectation, 
one can write ([F'] 2 ) = [</ r ')] 2 , so that equation (65) becomes 

MN) =2jdJ [<F'>] 2 . (66) 

Now the dependence of the term (F') with N follows from equa¬ 
tion (52), so that one can write 

hi(N) = ^ , (67) 

where a 2 is an amplitude independent of N. This scaling is a pre¬ 
diction from the Balescu-Lenard formalism. If the secular evolu¬ 
tion observed in S12 simulation was a Vlasov-only evolution, i.e. 
a collisionless evolution, one would expect a scaling of h 2 such 
that dh 2 /dN-0. 

One may now compare these predictions to the scalings ob¬ 
tained from A'-body runs. We considered number of particles 
given by Ne{8 , 12, 16, 24, 32, 48, 64) x 10 s , and for each of 
these values of N, we ran 32 different simulations with differ¬ 
ent initial conditions while using the A'-body code described 
in section 5. For each value of N, one may study the function 





. 8x10 5 
. 12x10 5 
• 16x10 s 
. 24x10 5 
. 32x10 5 
. 48x10 5 
. 64x10 5 


Fig. 5: Illustration of the behavior of the function t\->h(t,N) 
defined in equation ( 57 ), for an active fraction £ = 0.5, when 
averaged on 32 different realizations for particles numbers 
Ne {8, 12, 16, 24, 32, 48, 64)xl0 5 . To compute h(t,N), we binned 
the action-space domain (/<j, J r ) = [0;2.5]x[0 ;0.2] in 100x50 regions. 
The values of h(t,N) have also been uniformly renormalized so as to 
clarify this representation. The dots corresponds to the snapshots of 
the simulations for which h(t,N) was computed, whereas the lines 
correspond to second-order fits. As expected, the smaller the number of 
particles, the noisier the simulation and the larger h(t, N). 


Log(MN)) 



Log(MN)) 



Log(N) 


Fig. 6 : Top panel: Illustration of the behavior of the function 
log(/V) i-> log(/io(AO), where N has been rescaled by a factor 1(L 5 
so as to simplify the representation. The dots correspond to com¬ 
puted values thanks to figure 5 , while the line corresponds to a linear 
fit, which takes the form log(/i 0 (AO)~ 11.75-1.02 log(A0. The coeffi¬ 
cients / 70 (A) have been uniformly renormalized so as to clarify this 
representation. Bottom panel: Similar representation for the behavior 
of the function log(A)i->log(/! 2 (A)), whose linear fit takes the form 
log(/t 2 (A))= 12.36-1.91 log(A). Similarly, the coefficients h 2 (N) have 
been uniformly renormalized so as to clarify this representation. 


t\->h(t,N), as illustrated in figure 5. Once the behavior of the 
function 1 1 -» h(t, N) is known, one can fit to these parabolas as in 
equation (59), so as to determine the behavior of the functions 
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N I—> h{)(N) and li 2 {N). The dependence with N of these coeffi¬ 
cients is illustrated in figure 6. From the top panel of figure i, 
we recover the scaling of ho(N) derived in equation (62) due to 
Poisson shot noise. The bottom panel of figure 6 displays the 
scaling h 2 (N)ccN~ 191 . Given the finite number of simulations 
considered and the uncertainties in the fits, this is in good agree¬ 
ment with the result presented in equation (67). This scaling of 
Ji 2 (N) with N therefore confirms the relevance of the Balescu- 
Lenard formalism in describing the secular evolution of S12 sta¬ 
ble Mestel disc. Specifically, as explained below equation (67), if 
the features observed in S12 simulation had only been the result 
of a collisionless mechanism, one would not have observed such 
a scaling of /;i(A0 with N. This scaling confirms that the secular 
evolution of S12 stable Mestel disc is the result of a collisional 
evolution seeded by the discrete nature of the system and the ef¬ 
fect of amplified distant resonant encounters. Another probe of 
the collisional scaling, which allows to get rid of Poisson shot 
noise as present in equation (62), is described in Appendix F. 


5.3. Scaling with £ 


Since the novelty of the Balescu-Lenard formalism is to capture 
the effect of gravitational polarization, we now further compare 
qualitatively the prediction from section 4 with the results ob¬ 
tained from numerical simulations, by studying the impact of 
the active fraction £ of the disc on the observed properties of 
the secular diffusion. Indeed, as detailed in section 4.1 , the disc 
considered in S12 had an active fraction of £ = 0.5, so that only 
one half of the potential was due to the active component. If one 
increases the active fraction of the disc, one will increase the 
strength of the self-gravitating amplification, and therefore ac¬ 
celerate the secular evolution of the disc, while still remaining 
in a regime of collisional evolution. Therefore the scaling of I12 
with N given by equation (61 ) will remain the same, but the pref¬ 
actor c* 2 (£) will increase because the secular evolution will be 
amplified via a more efficient polarization. The dependence of 
<*2 with £ can be both measured from A'-body simulations but 
also predicted using the Balescu-Lenard formalism via the cal¬ 
culations presented in section 4.2. Let us consider the same sets 
of simulations as in section 5.2, so that the number of particles 
were given by Ne {8, 12, 16, 24, 32, 48, 64}x 10 5 , and for each 
of these values of N, 32 different simulations with £ = 0.6 were 
performed, in order to carry out ensemble averages. 

The equivalent of figures 5 and 6 for £ = 0.6 is illustrated in 
figure 7 . Even for £ = 0.6, one finds that the function t^ hit, N) 
follows a parabola given by equation (59). One also recovers 
the predicted scalings with N of the functions Nt-^hu(N) and 
(VhA 2 (W), respectively representing the initial Poisson shot 
noise of the sampling and the collisional scaling of the Balescu- 
Lenard secular evolution. As expected, when the active frac¬ 
tion of the disc is increased the secular evolution is fastened. 
Thanks to these fits, one can study the dependence of the ratio 
<* 2 (£ —0.6)/ar2(£ = 0.5), as defined in equation (67), both from nu¬ 
merical simulations as described in figures 6 and 7 and from the 
Balescu-Lenard equation using the matrix method described in 
section 4. 

From the fits of N^>!i2(N) from figures 6 and 7, one 
can write log(/z 2 (iV)) = 6.40-1.91 (log(A)-3.12) for £ = 0.5 and 
log(7z2(TV)) — 9.76 —1.84(log(7V)-3.12) for £ = 0.6, where we 
shifted the intercept of the fits to correspond to the center of 
the considered region log(7V)e [log(8); log(64)]. One therefore 
obtains the ratio 


<*2(0.6) 

<*2(0.5) 


= exp [9.76-6.40] = 29. 

NB 


( 68 ) 


ft(t,N) 



• 8x10 s 

• 12x10 s 

• 16x10 s 

• 24x10 s 

• 32x10 s 

• 48x10 s 

• 64x10 s 


log(/? 0 (N)) 



Log(N) 


Log(MN)) 



Log(N) 


Fig. 7: Top panel. Illustration of the behavior of the function 1 1 -» h(t, N) 
for an active fraction £ = 0.6, using the same conventions as in figure 5. 
As expected, for an increased active fraction, the secular evolution of 
the system is fastened leading to a faster increase of the distance h(t, N). 
Middle panel: Behavior of the function log((V)i->log(/7o((V)), for an ac¬ 
tive fraction £ = 0.6, using the same conventions as in figure 6. Its lin¬ 
ear fit takes the form log(/7o(/V)) =11.90- 1.07 log(lV). One recovers the 
expected scaling of the Poisson shot noise sampling derived in equa¬ 
tion (62). Bottom panel. Behavior of the function log(A0>->log(/t2(A0), 
for an active fraction £ = 0.6, using the same conventions as in figure 6. 
Its linear fit takes the form log(/i2(A0) = 15.48-1.84 log(iV). One recov¬ 
ers the expected collisional scaling with N obtained in equation (67). 


One may now compare this /V-body measurement to the 
same measurement performed via the Balescu-Lenard formal¬ 
ism. Following equation (66), one obtains that this ratio is given 
by 


»,<?,) 

Jd2 (diver*)] 2 


where yf ot stands for the secular diffusion flux at f = 0 1 with an 
active fraction £. The value of <* 2 (£ = 0.5) can be determined via 
figure t, while the secular diffusion flux determined for £ = 0.6 
is illustrated in figure 8. Thanks to the contours presented in fig- 
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Fig. 8: Map of VdivfTqot) with an active fraction £ = 0.6 and using the 
same conventions as in figure . The contours are spaced linearly be¬ 
tween the minimum and the maximum of fVdiv(!F tot ). The maximum 
value for the positive blue contours corresponds to Vdiv(!Ftot)-4200, 
while the minimum value for the negative red contours corresponds to 
MiivCFto,)^ -3200. As expected, one recovers that when the active frac¬ 
tion of the disc is increased, the susceptibility of the disc is increased, 
so that the norm of Vdiv(!Ft 0 t) gets larger and the secular diffusion is 
fastened. 


4 

A. 


0 0.5 1.0 1.5 
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ures 4 and 8, one can perform the same estimation as in equa¬ 
tion (68) starting from the Balescu-Lenard predictions. In or¬ 
der to focus on the contributions associated with the resonant 
ridge, the integrals on J in equation (69) were performed for 
./^ e [0.5; 1.2] ande [0.06 ; 0.15]. We measured 


homogeneous system. Such a collisional evolution is only rele¬ 
vant for stable systems, i.e. systems assumed to be stable in the 
Vlasov sense. Because it has been obtained via a Taylor expan¬ 
sion of the dynamics at the order 1 /V in the number of particles, 
it remains valid only for secular timescales of the order Vfo, with 
to the dynamical time. 

On such secular timescales, a Balescu-Lenard evolution can 
lead to two different outcomes. On the one hand, if the system 
remains stable during its entire evolution, the Balescu-Lenard 
equation will tend towards a 1 /V-stationary state . Once such a 
stationary state of evolution has been reached, the dynamics is 
then governed by the next order kinetic effects in 1 /V 2 , which 
are not captured by the Balescu-Lenard equation. On the other 
hand, the Balescu-Lenard collisional evolution may lead also to 
a destabilization of the system. Indeed, the long-term effects of 
the collisional diffusion, because they lead to an irreversible dif¬ 
fusion of the DF, may change its current state w.r.t. the colli¬ 
sionless (Vlasov) dynamics. After a slow and stable evolution 
sourced by collisional 1 /V effects, the system may then become 
unstable with respect to collisionless dynamics, which becomes 
the main driver of its later-time evolution, as was suggested by 
Sellwood (2012). S12 observes an out-of-equilibrium transition 
between the 1 /V Balescu-Lenard collisional evolution and the 
collisionless Vlasov evolution. 

One can illustrate such a transition using the V-body sim¬ 
ulations presented previously. In order to capture the change of 
regime of evolution within the disc (collisional vs. collisionless), 
for a given value of the number of particles, we define the quan¬ 
tity T. 2 (t,N) as 

^2(1, N) = ( J^d/e R d<b I star (t, V, R, <!>) e~ l2 ^ = (/£ e- lU j , (71) 


(* 2 ( 0 - 6 ) 

(*2(0.5) 


^ 42 . 

BL 


(70) 


Despite the difficulty of this measurement which required to 
consider a much more sensitive disc with £ = 0.6, the ratios of 
(*2 measured either via direct /V-body simulations as in equa¬ 
tion (68) or via application of the Balescu-Lenard formalism in 
equation (70) are within the same order of magnitude. As a con¬ 
sequence, one indeed checks that the Balescu-Lenard equation 
is able to correctly capture the relative effect of the disc suscepti¬ 
bility on the characteristics of the collisional secular diffusion. 
The strong consequence of modifying the active fraction, ob¬ 
served both in equations (68) and (70), illustrates the relevance 
of the self-gravitating amplification in determining the typical 
timescale of secular diffusion of the system. 


5.4. Late-time evolution 

The predictions of the Balescu-Lenard secular diffusion flux pre¬ 
sented in section 4 were only applied for the initial time of evolu¬ 
tion, i.e. for the estimation of T\ (]l (t-Q ¥ ). The V-body simula¬ 
tions presented in section 5 allowed us to verify the appropriate 
scaling of the response of the system with the number of par¬ 
ticles for the initial time of evolution as illustrated in figure 6. 
Using the Balescu-Lenard formalism to probe the late secular 
evolution of the system would require to evolve iteratively equa¬ 
tion (2) over secular times. Such iterations are clearly beyond the 
scope of this first paper, however the use of V-body simulations 
allows us to start probing now such late times of evolution. 

As discussed in section 2, the Balescu-Lenard equation de¬ 
scribes the long-term evolution of a discrete self-gravitating in¬ 


where as in equation ( 57 ), the operator (•) corresponds to the 
ensemble average, approximated here with the arithmetic aver¬ 
age over the p = 32 different realizations of simulations for the 
same number of particles V. The radii considered are restricted 
to the range Re[R{ n f; f? sup ] = [1.2 ; 5], where the active surface 
density of the disc is little affected by the inner and outer tapers. 
Finally, to obtain the second equality in equation (7 1), as in equa¬ 
tion (C.2), we replaced the active surface density of the disc by a 
discrete sum over all the particles of the system, where the sum 
on n is restricted to all the particle whose radius lies between R ul f 
and /? sup , while their azimuthal phase was written as ([>„. Such 
a quantity allows us to probe easily the presence of strong non- 
axisymmetric features within the disc. During the initial Balescu- 
Lenard collisional evolution of the system, one expects low val¬ 
ues of J. 2 - Indeed, during this evolution, one relies on the phase 
averaging approximation, which assumes that F = F(J, t), so that 
the DF of the system does not depend on the angles 9. During 
this collisional phase, 5L still remains non-zero because the sys¬ 
tem develops transient spiral waves, which sustain the secular 
evolution. On the long-term, this collisional evolution, through 
an irreversible diffusion of the DF, leads to a destabilization of 
the system. Eventually, the dynamical drivers of evolution are 
not any more discrete resonant collisionless effects but exponen¬ 
tially growing dynamical instabilities. In this regime of collision¬ 
less unstable evolution, one expects much larger values of 5L, be¬ 
cause of the appearance of strong non-axisymmetric bars within 
the disc. This bifurcation between these two regimes of diffu¬ 
sion is illustrated in figure 9, through the behavior of the func- 

2 Boltzmann DF of the form exp[-/3 H{J )], when physically reachable, 
are obvious stationary states of the Balescu-Lenard equation. 


Article number, page 12 of 21 










J. B. Fouvry et al.: Secular diffusion in discrete self-gravitating tepid discs 


« 1/2 Z 2 (t,N) 



- 8x10 s 

- 12x10 s 

- 16x10 s 

- 24x10 s 

- 32x10 s 

- 48x10 s 
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Fig. 9: Behavior of the function fi-> VAX 2 (f) as defined in equa¬ 
tion (71), for various values of the number of particles. The prefactor 
VA has been added so as to mask Poisson shot noise, allowing the ini¬ 
tial values of VA E 2 to be independent of A. It illustrates the bifurca¬ 
tion between the initial Balescu-Lenard collisional evolution, for which 
low values of E 2 are expected and the collisionless Vlasov evolution for 
which the system is no more axisymmetric leading to larger values of E 2 . 
As expected, the larger the number of particles, the later the transition. 


tion th*I, 2 (t,N) . One can similarly observe this transition di¬ 
rectly by looking at the active surface density X t (/7,0, t) for these 
two different regimes. This is illustrated in figure 10, where one 
recovers that in the late time collisionless regime of evolution, 
the galaxy becomes strongly non-axisymmetric. S12 found that 
just after the disc becomes unstable, the pattern of the spiral re¬ 
sponse is consistent with the ILR frequency corresponding to 
the ridge . Hence the phase transition observed in figure 10 is 
driven by all the free energy available in a cold disc, which via 
spiral transients secularly heats the disc, but only along a very 
tight resonant direction. This in turn leads the disc towards an 
orbital instability, transverse to the resonance (via the direct az¬ 
imuthal analogue to the two stream instability in plasma physics, 
Lynden-Bell 1979; Pichon 1994). Qualitatively, one expects that 
the more massive and the narrower the ridge, the larger the num¬ 
ber of orbits trapped in ILR resonance with little relative az¬ 
imuthal dispersion , and the earlier the instability (Penrose 1960; 
Pichon & Lynden-Bell 1993). 

In closing, it is quite striking that an isolated galactic disc, 
fully stable in the mean field sense, will, given time, drive it¬ 
self through two-point resonant correlations towards instability, 
demonstrating the extent to which such cold systems are truly 
secularly metastable. 

6. Conclusion 

Most astrophysical discs formed through dissipative processes 
and have typically evolved over many dynamical times. Even 
in isolation, the long range force of gravity allows their com- 

3 A similar dynamical phase transition has been observed (Campa et al. 
2008) in a toy model of systems with long-range interactions called 
the Hamiltonian Mean Field (HMF) model. During the slow collisional 
evolution, because of finite-A effects, the distribution function of the 
system changes with time. In certain cases, the system may become dy¬ 
namically (Vlasov) unstable and undergo a rapid phase transition from 
a homogeneous phase to an inhomogeneous phase. This phase transi¬ 
tion can be monitored by the magnetization (see Fig. 1 of Campa et al. 
(2008)) which is an order parameter playing a role similar to X 2 (t, A). 

4 One could also check that the disc’s distribution function corresponds 
at that stage to an unstable configuration, using the matrix method de¬ 
scribed in Appendix C. 


Fig. 10: Illustration of the active surface density E t for a A-body run 
with A = 8xl0 5 , restriced to the range R< 6. Top panel. Active surface 
density X, at an early time t = 60, for which the galaxy remains globally 
axisymmetric. In this regime, the dynamics of the system is collisional 
and governed by the Balescu-Lenard equation (2). Bottom panel. Ac¬ 
tive surface density E t at a much later time f = 2400. The galaxy is then 
strongly non-axisymmetric. In this regime, the dynamics of the system 
is collisionless and governed by the Vlasov equation. 

ponents to interact effectively through resonances, which 
given time may drive them secularly towards more likely 
equilibria. Such processes are captured by recent exten¬ 
sions of kinetic theories rewritten in angle-action variables 
(Heyvaerts 2010; Chavanis 2012). Solving these equations 
provide astronomers with a unique opportunity to quantify the 
induced secular angular momentum redistribution within these 
discs (Lynden-Bell & Kalnajs 1972) over cosmic timescales. 
While challenging, the numerical computation of the corre¬ 
sponding diffusion and drift coefficients is as demonstrated 
within reach of a relatively straightforward extension of the 
so-called matrix method (Kalnajs 1976), which computes the 
orbital response of self-gravitating discs using quadratures and 
linear algebra. 

Paper I presented asymptotic expressions in the tightly 
wound limit and provided a qualitative insight into the phys¬ 
ical processes at work during the secular diffusion of a self- 
gravitating discrete disc. Conversely, in this paper, we computed 
numerically the drift and diffusion coefficients of the inhomo¬ 
geneous Balescu-Lenard diffusion for such infinitely thin stellar 
discs. The self-gravity of the disc was taken into account via 
the matrix method, validated on unstable Mestel discs. We com¬ 
puted the divergence of the flux density in action space, div(lFtot)- 
Swing amplification was shown to provide a significant boost 
for the diffusion timescale, which now matches the numerically 
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measured one. These computations are the first exact calcula¬ 
tion of the Balescu-Lenard diffusion and drift coefficients in 
the context of inhomogeneous multi-periodic systems. They cap¬ 
ture the essence of self-induced evolution (nature), which should 
compete with environmentally induced evolution (nurture). We 
then compared these predictions to idealized numerical simu¬ 
lations of stable stationary and truncated Mestel discs sampled 
by pointwise particles, which were evolved for hundreds of dy¬ 
namical times. Using ensemble averages of our A'-body runs, 
we also identified a clear signature of the Balescu-Lenard pro¬ 
cess in the scaling of the diffusion features with N and <f, the 
fraction of the mass within the disc. As originally identified 
by Goldreich & Lynden-Bell (1965); Julian & Toomre (1966) in 
the context of their linear response, the susceptibility of cold 
self-gravitating discs plays a critical role for their secular evo¬ 
lution as it is squared in the Balescu-Lenard equation, which 
boosts considerably the effect of discreteness. Indeed, both the 
numerical experiments and our computation of the fluxes show 
that N e fi ~N/10 4 , which is consistent with the predicted rescal¬ 
ing in 1 /I) 2 (it was shown forty years ago that for a Mestel 
disc 1/2) ~ 10 2 , depending on the exact temperature of the disc, 
Toomre 1981). 

Jointly with paper I we now have a qualitative and quanti¬ 
tative understanding of the initial secular orbital diffusion pro¬ 
cess induced by the discreteness of galactic discs. Our qualita¬ 
tive understanding allows us to identify the role played by the 
square susceptibility in boosting the diffusion. Our quantitative 
agreement in both amplitude, position, width and scaling of the 
induced orbital signatures strongly suggests that secular evolu¬ 
tion is indeed driven by resonances as captured by the Balescu- 
Lenard formalism, and that it does not depend on the initial 
phases of the disc (since the matching Balescu-Lenard fluxes are 
phase averaged). It demonstrates that this equation initially re¬ 
produces the observed evolution of self-gravitating discs driven 
by resonant two-point correlations beyond the mean field approx¬ 
imation. 

The next step will be to evolve iteratively equation (2) over 
a Hubble time, and compare with the result of A'-body simula¬ 
tions. One should also model it jointly with an externally induced 
orbital diffusion (Fouvry et al. 2015b) arising from e.g. a (possi¬ 
bly anisotropic) cosmic environment (Codis et al. 2012, 2015) 
so as to assess which process dominates. At the technical level, 
the Balescu-Lenard formalism should be used to (in)validate 
A-body integrators accuracy over secular timescales. There are 
indeed very few analytical predictions on which to calibrate 
A'-body experiments in this regime. Such an exploration would 
also allow us to get a better grasp of the impact of the numeri¬ 
cal parameters used in the A'-body integration (such as timestep, 
mesh size or softening length) on the long-term dynamics of the 
system. 

Beyond the application described in this paper, the Balescu- 
Lenard formalism may in the future also be numerically im¬ 
plemented to describe for instance the secular diffusion of gi¬ 
ant molecular clouds in galactic discs (which in turn could play 
a role in migration-driven metallicity gradients and disc thick¬ 
ening), the secular migration of planetesimals in partially self- 
gravitating proto-planetary debris discs, or even the long-term 
evolution of population of stars and gas blobs near the Galactic 
centre. In 3D, assuming spherical symmetry, its implementation 
could be useful to describe spherical systems dominated by ra¬ 
dial orbits, or the secular evolution of tidal debris in our possibly 
flattened galactic halo using Stackel potentials. 
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Appendix A: Kalnajs basis 


Appendix B: Calculation of S 


We now detail the properties of the basis introduced in Kalnajs 
(1976) to describe 2D discs . The basis will depend on two pa¬ 
rameters: an index and a scale radius /- K a eR + • In all the 

upcoming formula of this section, in order to shorten the nota¬ 
tions, we will write r for the dimensionless quantity r/r Ka - As 
introduced previously in equation (16), the basis elements will 
depend on two indices: the azimuthal number £ and the radial 
index n. One should note that we have £ > 0 and n> 0. The radial 
component of the potential elements are then of the form 


Vg 


mr)=-^ w nkK a 

r Ka 


k n 

JVH Vi 


t, n) r y t y a Ka (kKa,{,n,iJ) r 
i =0 j =0 


2i+2j 


(A.l) 


We now detail how the analytical function N from equation (32) 
may be computed. In order to ease the implementation of its com¬ 
putation, we rewrite N in an undimensionnalized way as follows 

H(a g ,b g ,c g ,a h ,b h ,c h ,T],Ar) = 



The radial component of the density elements is given by 
£>i(r) S(k Ka , £, n) (1 -r 2 )^- 1/2 / 

^ Gr L 

n 

xJ]/5 Kd (k Ka J,nJ)(l-r) j . (A.2) 

j =o 

In equations (A ) and (A. 2), the coefficients P(k,£,n) and 
S(k, £, n) are defined by 


P(kJ,n) = 


\[2k+C+2n + (\l2)]Y[2k+£+n+ai2)] 

\ r[2Ar+/n-i] r 2 [^-ni]r[H+i] 

1 1/2 

xf[f+fl+(l/2)] 


(A.3) 


and 


S(k,£, n) = 


T[k+1] 


zr r[2A:+1] T[^+(l/2)] 
f [2k+£+2n+(l/2)\ Y[2k+n + \] T[2k+£+n- 


( 1 / 2 )]) 


1/2 


r[^+«+(i/2)]r[n+i] 


(A.4) 


Finally, in equations (A3 ) and (A.2), we have also introduced 

n . . . [—&]; [^+(1/2)]/ [2k+£+n+(\/2)\j 

avJk, c, n, i , /) =- 

x[i+£+(V2)]j[-n\j, (A.5) 


and 


Pkz( k, £, n,j) 


[2k+£+n+(l/2)]j[k+l]j[-n]j 

[2*+l];[* + (l/2)];[l]; 


(A.6) 


In the two previous expressions, we introduced the rising 
Pochhammer symbol [a], defined as 


M; = 


1 if / = 0, 

a (a+ 1) ... (fl+n-1) if i> 0. 


(A.7) 


5 See also Earn & Sellwood (1995 j for a similar rewriting of the basis 
normalizations. 


where we assumed that a g ,ai ,^0 and used the change of vari¬ 
ables x = x p / A r and y = x a / A r. We also defined the dimensionless 
function Nd as 


N d (&, c,e,f,rf) = 


1 J- 1 

2 2 


d.vdzy 


1 +bx+cy 
1 +ex+ fij+irj 


(B.2) 


To compute this integral, we may now exhibit a function G(x, y) 
such that 


d 2 G 


1 +bx+cy 


(B.3) 


dxdy l+ex+fy+ir) 

One possible choice for G is given by 

G(x, y) = -^^\og[e 2 x 2 +2e(fxy+x)+f~y 2 +2fy+y 2 +1] 

\y>f(e 2 x 2 -(fy+irj+ l) 2 )+2e f(ex+ir/+ \ )-ce(ex+h]+ 1) 2 | 

i r 7r ,\ex+fy+\V\ 

2 e^f 2 12 rj I 

x | bf(e 2 x 2 - (fy + irj +1 ) 2 )+ 2ef(ex+ir\ +1) - ce(ex+ irj +1 ) 2 1 
+ |/He +b(2ex+fij+2ir]+2)) 

+ce(2ex- fy+2irj+2)+2ef(cy+2) log[ex+fy+ir/+ 1] | (B.4) 


One should note in the previous expression the presence of a 
complex logarithm and a tan -1 . However, because e, /, rj e R and 
/ 74 O, one can easily show that the arguments of both of these 
functions never cross the usual branch-cut of these functions 
{lm(z) = 0; Re(z)<0). As a consequence, the expression ( .2) 
can immediately be computed as 

Nz) = G[i,i]-G[2,-i]-G[-i,i]+G[-i,-2]. (B.5) 


Appendix C: Response Matrix and N -body 
validations 

The computation of the response matrix as described in sec¬ 
tion 3 was validated by recovering the results of the pioneer 
work of Zang (1976), extended in Evans & Read (1998a,b), and 
recovered numerically in Sellwood & Evans (2001). These pa¬ 
pers predicted the precession rate <jJo = m < pQ. p and growth rate 
y-s of the unstable modes of a truncated Mestel disc similar 
to the stable one described in section . To build up an unsta¬ 
ble disc similar to the ones considered in these previous works. 
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one has to consider a fully active disc, so that £ = 1. So as to have 
Q -1 (Toomre 1964), the velocity dispersion within the disc will 
be given by q — 6 , where the parameter q has been introduced in 
equation ( 8 ). Finally, a last parameter one can tune in order to 
modify the properties of the disc is the truncation index of the in¬ 
ner tapering v t defined in equation (50). While looking only for 
m,p = 2 modes, we considered three different truncations indices 
given by v t = 4, 6 , 8 . To compute the response matrix, we used 
the same numerical parameters as described in section 1.2. Look¬ 
ing for unstable modes amounts to looking for complex frequen¬ 
cies u) = u>o + irj, such that the response matrix M(wo ,']) from 
equation (24) possesses an eigenvalue equal to 1. Such a complex 
frequency is then associated with an unstable mode of pattern 
speed u)q and growth rate q. To determine the growth rate and 
pattern speed of the unstable modes, we relied on Nyquist con¬ 
tours similarly to the technique presented in Pichon & Cannon 
(1997). For a fixed value ofone can study the continuous com¬ 
plex curve wot—>det [I—M(wo, 77 )]. Because for 77 —»+ 00 , one has 
|M(«o, 77 H —> 0, the number of windings of this curve around the 
origin gives a lower bound on the number of unstable modes with 
a growth rate superior to 77 . By decreasing the value of 77 , one can 
then determine the largest value of 77 admitting an unstable mode, 
and therefore the most unstable mode of the disc. The Nyquist 
contours obtained for the truncation index v t = 6 are illustrated in 
figure C. 1 , while the measurements are gathered in figure C.5. 


lm[det[l-M(w)]] 

0.02 r 


0.01 



'-0t)L/ / ! 
/ /- 0 / 0 V 

/lit 

1 ' ' L 
' ' -‘O.Cft 


\ 0.01 


“ 2 Re[det[l-M(w)]] 


Log |det[l-M(w)]| 



0.24 

0.22 
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Fig. C. 1: Top panel: Zoomed Nyquist contours in the complex plane of 
ojq 1 —> det [I— M(tu 0 , 77 )] obtained via the matrix method for a truncated 
Mestel disc with v, = 6 and q- 6, looking for m,p = 2 modes. Each contour 
corresponds to a fixed value of 77. For a growth rate of 77 = 0 . 20 , one 
can note that the contour crosses the origin, which corresponds to the 
presence of an unstable mode. Bottom panel. Illustration of the function 
oio i-> log | det [I—M(oio, 77)] | for the same truncated Mestel disc. Each 
line corresponds to a fixed value of 77. Such a representation allows to 
determine the pattern speed oi 0 = 777^0^ = 0.94 of the unstable mode. 


Once the characteristics (a>o, 77 ) of the unstable modes have 
been determined, one can study in the physical space the shape 
of the mode. Indeed, for to-aio+iq, one can compute M(o»o, q), 
and numerically diagonalize this matrix. One then considers its 
eigenvector Z m 01 i e (of size 77 max , where 77 max is the number of 
basis elements considered) associated with the eigenvalue almost 


equal to 1. The shape of the mode is then immediately given by 


Smode(W <t>) = Re 


de^W) 


(C.l) 


where 'L ip] are the considered surface density basis elements. The 
shape of the recovered unstable mode for the truncated v t = 4 
Mestel disc is illustrated in figure C.2. 



Fig. C.2: Unstable mode for the truncated v t =4 Mestel disc recovered 
via the matrix method presented in section 3. Only positive contour lev¬ 
els are shown and are spaced linearly between 10% and 90% of the 
maximum norm. The radii associated with the resonances ILR. COR 
and OLR have been represented, as given by aJo = m-Sl(R m ), where 
the intrinsic frequencies Sl(R) = (Q. l/ ,(R), k(R)) have to be computed 
within the epicyclic approximation. For a Mestel disc, they are given 
by Q^R) = T 0 /R and k(R) = V2 Cl^R). 


The same unstable modes were also used to validate the 
TV-body code presented in section 5. To run these simulations, 
we used the same samping technique as described in section E. 
In order not to be significantly impacted by the absence of a quiet 
start sampling (Sellwood 1983), for each value of v t , the mea¬ 
surements were performed with simulations of 20 M particles. As 
observed in Sellwood & Evans (2001), the appropriate setting of 
the parameters of the TV-body code are crucial to recover cor¬ 
rectly the unstable modes of a disc. We considered a grid made 
with TV mes h = 120 grid cells, while using a softening length equal 
to s-Ri/60. As described, in section 5, we similarly restricted 
the perturbing forces only to the harmonic sector 771 ^ = 2, using 
TV ri ng =2400 radial rings, with = 720 azimuthal points. In or¬ 
der to extract the properties of the mode present within the disc, 
one may proceed as follows. For each simulation snapshot, one 
can estimate the active surface density within the disc via 

£star(*, 0 = f X > (C.2) 

i 

where the sum on i is made on all the particles of the simula¬ 
tion and Xj(t) is the position of the 7 th particle at time t. Such a 
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surface density can be decomposed on the basis elements from 
equation (3), under the form 

^starC^? 0 — z (C.3) 

p 


where the sum on p is made on all the basis elements I. tp> con¬ 
sidered. The effective basis elements used during our measure¬ 
ments are the same as the ones used in the matrix method from 
section 1.2. Thanks to the biorthogonality property from equa¬ 
tion (3), the coefficients b p (t ) can be immediately determined as 

b p (t) = ~fdx Wjc, t) ^ p)t (x) = -p J] ^(*.-(0) • (C.4) 

As we are looking for unstable modes within the disc, we expect 
to have b p (t) ocexp[-i(u>o + ir/)t], where ojq = m$£l p is the pattern 
speed of the mode and j- its growth rate. As a consequence, one 
immediately obtains that 


d Re(log(7? p (f))) 

dr 


d Im(log(fc p (f))) 
dr 


-wo, 


(C.5) 


if one is sufficently careful with the branch-cut of the complex 
logarithm. Such a linear scaling with t of Re(log(/?p(f))) and 
Im(log(Z? p (r))) is therefeore the appropriate measurement proce¬ 
dure to use in order to estimate the growth rate and pattern speed 
of the unstable modes of these truncated Mestel discs. These 
measurements for the various values of the truncation index v t 
are illustrated in figure C.3. Once the basis coefficients b p (t) 


Re(Log(fc 0 (t))) 



lm(Log(i> 0 (t))) 



Fig. C.3: Measurement of the growth rate q and pattern speed aio of 
the = 2 unstable mode truncated Mestel discs with a random velocity 
given by q = 6 for various values of the truncation index v t = 6, 8. The ba¬ 
sis coefficient plotted is associated with the basis element (£, n) = (2,0), 
using the same basis elements as for the matrix method in section 

have been determined, one can study the shape of the recovered 


unstable modes in the physical space. Indeed, similarly to equa¬ 
tion (C. ), the shape of the modes is given by 


^-mode(fC tf), 0 ~ Fc 




(C.6) 


In analogy with figure C.2, for which the unstable modes have 
been obtained via the matrix method, figure C.4 illustrates the 
unstable mode of the same truncated v t = 4 Mestel disc. 



Fig. C.4: Unstable mode for the truncated v t =4 Mestel disc recovered 
via direct A-body simulations as presented in section 5. Only positive 
contour levels are shown and are spaced linearly between 20% and 80% 
of the maximum norm. Similarly to figure C.2, the radii associated with 
the resonances ILR, COR and OLR have been represented. 


As a conclusion, the growth rates and pattern speeds obtained 
either via the matrix method or direct A-body simulations are 
gathered in table C.5. As observed in Sellwood & Evans (2001), 


Unstable m^-2 modes of truncated Mestel discs, q 

= 6. 


v t 

= 4 

Vt 

= 6 

v t = 

= 8 

Method 

Wo 

'7 

Wo 

t) 

w 0 

t] 

Linear Theory 

0.88 

0.13 

0.90 

0.22 

0.92 

0.27 

Matrix Method 

0.93 

0.11 

0.94 

0.20 

0.95 

0.24 

A-body 

0.99 

0.13 

0.79 

0.19 

0.89 

0.26 


Fig. C.5: Measurements of the pattern speed u> 0 = m^Q. p and growth 
rate q=s for unstable m^ = 2 modes of tapered Mestel discs. The ve¬ 
locity dispersion within these discs is characterised by q- 6, and the 
inner truncation power indices are given by v t = 4, 6, 8. The theoretical 
values were obtained from linear theory in Evans & Read (1998b). Our 
measurements were either performed via the response matrix method as 
in equation (24), or via direct A-body simulations, using the A-body 
integrator described in section 5. 


the recovery of the unstable modes characteristics from direct 
A-body simulations when performed for truncated Mestel discs 
is a difficult task, for which convergence to the values predicted 
through linear theory may be difficult. 
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Appendix D: Why swing matters? 

Let us investigate here briefly the importance of self-gravitation 
and the completeness of the projection basis in capturing the role 
of swing amplification. 


Appendix D.l: Turning off the self-gravitating amplification 

In order to investigate the role of the self-gravitating amplifica¬ 
tion, one may perform the same estimation as presented in fig¬ 
ure 4, while neglecting collective effects. When neglecting col¬ 
lective effects, i.e. when assuming that MsO, one recovers the 
inhomogeneous Landau equation (Chavanis 2013) which reads 


dF_ 

dt 


n(2n) d p 


dj\ 


//I] ^ 

)(" 


d_ 

dj\ 


dJi I A mx ,ltl2 (JuJiT 

d 


x 6 o(m v Tli-m 2 £l 2 )[mi- -m 2 - JJ^j F (J0 


• (D.l) 


Equation (D.l) involves the bare susceptibility coefficients 
\A,„i.m 2 (Ji* Ji)\ 2 , which can be equivalently defined (see Ap¬ 
pendix B of paper I) by 

A mu m 2 UuJT) = -£<( 7 ,)^) 

P 

= -Jd0id0 2 u(\x(0i,Ji)-x(0 2 , J 2 )\) e Km, 0i ~ m2e2) , (D.2) 


where u(x ) is the binary potential of interaction potential given 
by u(x)--G/\x\ for gravity. This estimation therefore does not 
require to estimate the response matrix from equation (24), but 
one still has to perform integrations along the resonant lines as in 
equation (44). Let us provide a first numerical implementation of 
this equation in the context of galactic dynamics. The contours 
of /Vdiv(7 r ^ [ j c ) are illustrated in figure D.L Comparing the maps 



0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 


J(p 

Fig. D.l: Map of AdivfFj 1 ”'), corresponding to the bare secular diffu¬ 
sion flux, using the same conventions as in figure . The contours are 
spaced linearly between the minimum and the maximum of AdivfiF^). 
The maximum value for the positive blue contours corresponds to 
A , div(? r |’o t re )-0.30, while the minimum value for the negative red con¬ 
tours is associated with IVdivCFjJjj” 5 ) - -0.50. This figure is qualitatively 
similar to the one obtained in figure 9 of paper I. 


of the dressed diffusion flux /Vdiv(7 - i„i) from figure and the 
bare diffusion flux AdivCF^ 6 ), allows to assess the strength of 


the self-gravitating amplification. As expected, when turning off 
the self-gravity of the system, one reduces significantly the sus¬ 
ceptibility of the system and therefore slows down its secular 
evolution, by a factor of about 1000. One may also remark that 
while the secular appearance of a resonant ridge in the dressed 
diffusion from figure 4 was obvious, the shape of the contours 
obtained in the bare figure D. do not emphasize as clearly the 
appearance of such a narrow resonant ridge. One can still remark 
that the structure of the bare contours obtained in figure D.l is 
similar to what was obtained in figure 9 of paper I, through the 
WKB limit of the Balescu-Lenard equation. One can finally note 
that the amplitudes of the bare divergence contours obtained pre¬ 
viously are similar to the WKB values obtained in paper I. As 
a consequence, the comparison of figures and D. 1 emphasizes 
that the strong self-gravitating amplification of loosely wound 
perturbations is indeed responsible for the appearance of a nar¬ 
row ridge, while also ensuring that this appearance is sufficiently 
rapid, as observed in the diffusion timescales comparison from 
equation (54). 

Appendix D.2: Turning off loosely-wound contributions 

As emphasized in the Introduction, the WKB limit of the 
Balescu-Lenard equation presented in Fouvry et al. (2015a) was 
not able to capture the mechanism of swing amplification, which 
involves unwinding perturbations. By considering a complete 
and global basis as in equation (16), we have shown in figure 4 
how the missing amplification from Fouvry et al. (2015a) could 
be recovered. Using the numerical method of estimation of the 
secular diffusion flux as presented in section 3, one can try to 
recover the results obtained within the WKB formalism by care¬ 
fully choosing the considered basis elements generically intro¬ 
duced in equation ( 16) and chosen to be by Kalnajs basis ele¬ 
ments as detailed in Appendix A. We recall that each basis el¬ 
ement depends on two indices: an azimuthal index £ and a ra¬ 
dial one n. Because in S12’s simulation perturbations were re¬ 
stricted to the harmonic sector = 2, one only has to consider 
basis elements associated with £= 2. Moreover, as illustrated in 
figure D.2, the larger n the radial index, the faster the radial vari- 



Fig. D.2: Illustration of the radial basis elements of the 1-2 Kalnajs 
basis elements for ^Ka -7, defined in Appendix A, which were used in 
the estimation of the Balescu-Lenard diffusion flux in section . As the 
radial index n increases, the basis elements get more and more wound. 


ation of the basis elements and therefore the more tightly wound 
the basis elements. So as to get rid of the loosely-wound basis 
elements which are the ones which can get swing-amplified, we 
perform a truncation of the radial indices considered. Therefore, 
we define the secular diffusion flux AdivfT^™) computed in 
the same way than /Vdiv(W t( , t ) as presented in section 4.2, except 
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that the basis elements are such that n cm <n<n m!a , with n cut — 2 
and H ma x = 8. By keeping only the tightly wound basis elements, 
one can therefore consider the same contribution as the one con¬ 
sidered in the WKB limit presented in paper I. The contours of 
tVdiv(!F^ t KB ) are illustrated in figure D.3. One can note that the 



0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 

J(p 


Fig. D.3: Map of A^div(7 r ^ t KB ), corresponding to the dressed secular 
diffusion flux, using the same conventions as in figure . In order to 
limit ourselves only to tightly wound contributions, the basis elements 
associated with the radial basis index ne {0,1 ) have not been taken into 
account. The contours are spaced linearly between the minimum and the 
maximum of iVdivkF^™). The maximum value for the positive blue 
contours corresponds to /Vdiv(Tj ot )-0.7, while the minimum for the 
negative red contours is associated with Ndivi'F™™) — -4.5. This figure 
is to be compared to figure 9 of paper I. 


values of the contours obtained in the map of /Vdiv(7 7 ^' 1 KI ’) illus¬ 
trated in figure D.3 are in the same order of magnitude as the 
ones which were presented in figure 9 of paper I in the WKB 
limit. The presence of positive blue contours of /Vdiv(‘F^ l KB ) is 
also in agreement with a secular heating of the disc (i.e. an in¬ 
crease of J r ). However, these contours do not display a narrow 
resonant ridge as was observed in S12 simulation or in figure 4. 


Appendix E: Sampling of the DF 

In order to use the /V-body integrator described in section 5, one 
has to sample the particles according to the DF given by equa¬ 
tion (51 ). We introduce the probability distribution function F sp , 
normalized to 1 and thanks to which the sampling is performed. 
This probability DF F sp is directly proportional to the active dis¬ 
tribution function F stal from equation (5 ), so that we may write 

F sp (E, L) = C sp L q exp[ -E/o*\ T inner (L) T outel (L ), (E. 1) 

where C sp is a normalization constant which will be deter¬ 
mined in the upcoming calculations. Because the mapping 
(E,L) i-> (./,•, from equation (12) is not a trivial one, we 
will not perform the sampling of the stars in the action space 
(J r , J,!,), but rather in the (E, L)-space. Moreover, one should 
pay attention to the fact that the DF F sp from equation (E. 1) 
is a probability distribution function in the (jc, v) -space, so that 
d 2 jtd 2 «.F sp (x, v) is proportional to the number of particles in the 
infinitesimal volume d 2 xd 2 u around the position (x, v). As we 
want to sample the particles in the (E, L)— space, we introduce 
the function h sp (E, L) such that dEdLh sp (E , L) is proportional to 


the number of particles in the volume dEdL around the location 
( E , L). One can now determine h sp (E, L ) as a function of F(E, L). 
Indeed, we have 


h sp (E 




= dxdv dv(E'-E) 6 D (L'-L) F S JE , L) 


= 2n 
— 2n 


-Jdr r j :h; r d(’ t 6 D (E'-E) 6 D (L'-L ) F sp (E, L ) 
rJdrJ'dv t S D (E'-E)F sp (E,L'), (E.2) 


using the fact that the tangential velocity satisfies u t = L/r. The 
last step is then to perform the change of variable v r —» E. One has 
v^=2(E-il/ M (r))-L 2 /r 2 , so that dv r = dE/yj2(E-i// M (r))-L 2 /r 2 . 
Because the radial velocity can be both positive and negative, 
equation (E.2) takes the form 


h sp (E\ U) = 4tt 


-M dE 


6 D (E'-E)F sp (E,U) 
V2 (E-if/ M (r))-L' 2 /r 2 


4 n~ 


Qi (E',Z/) 


-F sp (E',L'), 


(E.3) 


where we used the definition ( 3) of the radial intrinsic frequency 
Qi. One can then correctly normalize the probability distribu¬ 
tion /i sp and determine the value of the constant C sp from equa¬ 
tion (E. ). One should pay attention to the fact that in addition to 
the tapering functions Tinner and 7’ (mtcr from equations (50), we 
also assume that no stars have orbits that extend beyond R max . 
As a consequence, the allowed region in the ( E , L)-space has to 
satisfy two constraints. First of all, the angular momentum L stal 
has to satisfy 


Emin - 0 < T s tar ^ 7? max V() — L n 


(E.4) 


Then, for a given value of L star , one can show that the energy of 
the star E star must satisfy the constraint 


^min (^star) — *Am 


Vo 


+ 2 °<£ star < 


-^star 


= £ , max (T s tar)- (E.5) 


These two constraints allow to completely characterize the 
( E , L)-space on which the sampling (E, L) will have to be per¬ 
formed. One finally has to satisfy the constraint 


(z.) 

dL d E 

Enin ^EmmiL) 


dEh s JE,L) = 1 


(E.6) 


Given the parameters presented after equation (51), one can nu¬ 
merically determine the value of the constant C sp which reads 


C sp ^ 1.4723 x 1(T 15 . (E.7) 

We may now proceed to the sampling of the coordinates of the 
particles. Up to the sign of its radial velocity, one star is char¬ 
acterized by the set {E stas: , L stal ,R stal , (/> star }. Given that the initial 
state is axisymmetric, the azimuthal angle of the star can be uni¬ 
formly sampled between 0 and 2n. The next step is then to suc¬ 
cessively sample (L star , E star ) and finally /? star , using successive 
rejection samplings as we will now detail. 

The heart of the rejection sampling is as follows. Let us as¬ 
sume that we want to generate sampling values from a function 
fix), from which it is difficult to sample. However, we assume 
that we have at our disposal another distribution function g(x) 
from which the sampling is simple, and such that there exists a 
bound M > 1 satisfying fix) < M gix). The smaller M, the more 
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efficient the sampling. One then has to proceed as follows: sam¬ 
ple both a proposition x from g and u uniformly between [0; 1], 
One then applies the selection 

a _ /(•*) . (u<a ==> x is kept. ^ 

Mg(x) \u > a => x is rejected. 


so that one has fn oc ly. However, one should note that for r —> r p / a , 
one has /«(r)—>+ 00 , so that the rejection sampling cannot be 
used without considering a probability DF g which also di¬ 
verges for r-*r p / a . In order to get rid of these divergences, 
instead of sampling the variable r, we will sample the angle 
u G [—7r/2 ; tt/ 2], where we have defined the mapping n-> u{r) as 


In order to have an efficient sampling, one should try to consider 
a function g close to /. 

We may now directly sample (L, L ) thanks to this algorithm. 
The true sampling function is /(£,!,)= L sp from equation (E.2). 
The simple sampling function is g(E,L) LX 1, defined on the domain 
characterized by the constraints from equation ( ) and ( .5). 

When performing a rejection sampling with such an uniform 
gtE.i,)- in order to determine the bound M ( e,l), one only has to 
determine an uniform bound for f(E,E)- With the numerical val¬ 
ues introduced after equation (51), one can check that ftE.i.) is 
such that 


f(E,L)(E,L ) <1.4. 


(E.9) 


The final element required to be able to perform the rejection 
sampling with f(E,L) is to be able to draw uniformly candidate 
(E,L) in the domains defined by the constraints from equa¬ 
tion ( ) and (E.5), which is equivalent as sampling candidates 

(E,L) from the uniform probability distribution function gtEj.)- 
To perform this uniform sampling, since the constraints from 
equation (E.5) are expressed for a given value of L star , it is more 
natural to first draw L sta i- and then L star . The probability distri¬ 
bution according to which L star has to be drawn is of the form 
/z,oc(£ , max (L)-£ , min (L)). When correctly normalized, it reads 


ML) 


1 


2 L max Vo 


2 R 2 V 2 

v max y 0 


■log 


Lmax F() 


(E.10) 


To sample L from /}, we will use another rejection sampling by 
introducing the additional simple probability distribution func¬ 
tion gi defined as 

9l(L) = log 

ft is straightforward to check that /l<(3/2) gi, so that we may 
use the bound M L = 3/2 to perform the rejection sampling of L star . 
The final remark is to note that sampling L from gi is simple 
since its cumulative distribution function Gi = I d L'giiL') can 
be inverted so as to read 


kmnv Vi) _ 


(E. 11) 


G l \u) = fimaxVoeXp 


1 + W_1 



(E. 12) 


where W_i is the lower branch of the Lambert function W(x) for 
xg [-1/e; 0], With all these elements, the rejection sampling of 
L s tai- following fi from equation ( .10) can be performed. 

Once L star has been drawn, it only remains to sample uni¬ 
formly L s tar on the interval LM a r G [Lm m (L s t a i): L mLlx (L slar )], as 
given by equation (E.4). Thanks to these uniformly drawn can¬ 
didates (Lstar, Lstar) and the uniform bound from equation (E.9), 
one can perform the rejection sampling from the probability dis¬ 
tribution ftE,). 

For L s tar and L st ar succesfully sampled, one may then sample 
the radius /Liar using a similar rejection sampling. The radius has 
to be sampled according to the probability distribution /« given 
by 


Mr) = 


O \/n 

sj2(E-ifr M (r))-L 2 /r 2 ’ 


(E.13) 


r p +r a r a -r p 

r{u) = JL ^- + ^ JL ^(«), (E.14) 

so that one naturally has r(-7t/2)-r p and r(n/2)-r a . The prob¬ 
ability distribution function from which u has to be sampled is 
immediately given by 

Mu) = —Y~ cos (u)p r (r(u )). (E.15) 

Using the fact that the maximum of f, is reached for it-n 12, 
one can then sample u from f, using a rejection sampling with 
a uniform control probability distribution function g„(u) = l/7t. 
Once u is known, it only remains to compute L sta r = KMstar), so 
that the sampling of all the required quantities for one star has 
been performed. 

The final step of the sampling of the particles is to determine 
the physical coordinates of the particles ( x,v) associated with 
the set {L star , L st ar,L sta r, </> s tar}- These physical coordinates are the 
ones which will be given to the A'-body integrator. We draw uni¬ 
formly the sign of the radial velocity e r G{—1, 1). Because we are 
considering a disc made only of prograde stars, one immediately 
obtains that the radial and tangential velocities ty and are given 
by 

| y r = £ r ^2(L - (Lstar)) ~L” ta r/^part ’ (g 

ITt = L s tar /Lstar • 

The final step of the transformation to the ( x , n)-coordinates is 
then straightforward, since one naturally has 


X — Lgtar COS((/star) , 

U = Lstar Sin(</>star) , 

x ■ ,, , (E.17) 

V x — U r COS( 0 s tar) U| SIU( (/star) , 

Vy — U r sin(^star) + UtCOS(</. s tar)- 

One should note that the sampling procedure described pre¬ 
viously does not correspond to a quiet start procedure (Sellwood 
1983), which would allow a reduction of the initial shot noise 
within the disc, as briefly discussed in section C, with regard to 
the validation of the /V-body code. 


Appendix F: Another test of the scaling with N 

One difficulty with the measurement of the scaling with N pre¬ 
sented in section 5.2 is that one has to disentangle the con¬ 
tributions from the initial sampling Poisson shot noise present 
through ho from equation (62) and the effects due collisional 
Balescu-Lenard diffusion scaling through hj from equation (62). 
Indeed, Poisson shot noise leads to fluctuations of the system DF 
about its mean value. In order not to be sensitive to such fluctua¬ 
tions, one could only consider fluctuations sufficiently large, i.e. 
fluctuations caused by an effective secular diffusion rather than 
caused by inevitable Poisson fluctuations. As a consequence, by 
restricting ourselves only to large fluctuations, we can get rid of 
Poisson’s effects. We therefore define the function V(t,N) as 

V(t,N) = j'dJx[(F(t,J,N))-(F(t=0,J,N))<Cy\ , (F.l) 
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where we introduced a threshold Cy < 0. Here x [ ■ ] is a charac¬ 
teristic function equal to 1 if J,N))-(F(t=0, J , A))<Cy), 
and 0 otherwise. As a consequence, V(t,N) measures the vol¬ 
ume in action space of the regions (depleted from particles, since 
Cy<0) for which the mean DF has changed by more than Cy. 
For a sufficiently large value of the threshold Cy, such a construc¬ 
tion allows not to be polluted by Poisson sampling shot noise. 
For the initial times, as in equation (52), it is straightforward to 
study the scaling of V(t, N ) with t and N. Indeed, one can write 

m, J, N))-m = 0, J , A0> - At divCFtot) - ^ div((T tot ). (F.2) 

Introducing Vo- Jd/^[div(!Ftot)<Cy], one can rewrite equa¬ 
tion (F.2) under the form 

V(t,N) = jjVo- (F.3) 

Therefore, for a fixed value of N, one expects to observe a lin¬ 
ear time dependence of the function V(t,N ), as illustrated in fig¬ 
ure . In order to test the scaling of equation (F.3) with N, one 
may proceed as follows. Introducing a threshold value V t hoicU for 
each value of N, we define the associated threshold time f t hoid(A0 
as 


V(fthoidW),A0 = Vthoid. 


(F.4) 


Thanks to the scalings from equation (F.3), one immediately ob¬ 
tains that 


hhold(A0 - N 


Vthold 

Vo 


(F.5) 


Such a linear scaling of /thoid(A0 with A is a prediction from the 
Balescu-Lenard formalism and is nicely recovered in figure F.l. 


Appendix G: Distributed code description 

For the sake of reproducibility, which has been lacking in the 
context of the linear response of stellar systems, we distribute 
the linear matrix response code we wrote for this paper both as a 
Mathematica package (http: //www. iap. fr /user s/pi chon/ 
matrix-method/code/matrix-method.m), and a notebook 
(http://www.iap.fr/users/pichon/matrix-method/ 
code/matrix-method.nb). The functions therein allow for: 

• the determination as a function of ( r p , r a ) of the orbits quan- 
tites: E, L. J r , Q| and IT- 

• the construction of the 2D basis from Kalnajs (1976). 

• the computation of the Fourier transform w.r.t. the angles, i.e. 
the computation of < 3V^ 1 m2 „ P (/) from equation (21). 

• the calculation of the 2D response matrix via equation (33). 

It has been tested for the isochrone and the Mestel disc. 
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Fig. F.l: Top panel. Illustration of the behavior of the function 
t V(t,N ) from equation ( 1), when averaged on 32 different realiza¬ 
tions for particles numbers Ne{8 , 12, 16, 24, 32, 48, 64)xl0 5 , along 
with the associated linear fits. To compute V(t,N), we used the same 
binning of the action-space (7^, J r ) as in figure 5. As obtained in 
equation (F.3), one recovers that for a fixed value of N, the function 
1 1 —> V(t, N) is linear. The horizontal dashed line illustrates the thresh¬ 
old value {/hold f°r which the threshold time f thold is determined. Bottom 
panel: Illustration of the behavior of the function Aha r thold (A). As de¬ 
rived in equation (F.5), one recovers a linear dependence of Tthoid with 
N. 
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